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Lf. INTRODUCTION 


The current design trend toward higher tip speed, 
reduced engine weight, and higher pressure ratios has made 
the modern high speed fan and Compressor more susceptible to 
supersonic unstalled flacter. Unstalled PLutver, as 
distinguished from stall flutter, occurs at small angle of 
attack and with attached flows. In addition, it occurs at or 
near the design condition. Flutter itself is a self-excited, 
unsteady phenomenon which can occur when the exciting forces 
just egual the elastic forces of the turbomachinery blades. 
For a review of the problem of fan and compressor flutter, 


refer to Ref. 6. 


The original work in flutter analysis evolved from the 
flutter of aircraft wings, which was observed as early as 
1914. Theodorsen [Ref. 22] developed the most basic analytic 
treatment of unsteady, subsonic aerodynamics for an airfoil 
oscillating in pitch and plunge. This analysis dealt with 
what has become known as classical flutter, and showed that 
flutter was only possible when the natural frequencies of 
flexure and torsion coalesced as the relative velocity over 
the wing was increased. Each mode of oscillation by itself 


(pitch or plunge) was: stable. 


Cascade flutter differs from classical flutter in that 
the inertia and stiffness properties of a turbomachine blade 
differ greatly from those of a wing. AS a result, the 
critical flutter speed, the speed at which the onset of 
flutter occurs, is very much greater. Garrick and Rubinow 
mhet. 9 j extended the work of Theodorsen to include 
supersonic flow, and their analysis predicted that single 
degree of freedom torsional flutter was possible, although 


the bending node was still stable. Lane [Re.. 14] showed 
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that cascade flutter could be analyzed by considering only a 


Single oscillating blade. 


Work on cascades with subsonic leading edge locus is 
fairly recent. The first attempt was made by Gorelov 
[Ref. 10]. Verdon [{Ref. 23] obtained some of the first 
numerical results using a finite difference procedure. Brix 
and Platzer [Refs. 3,4] were also able to obtain results 
using a method of characteristics approach. Another 
solution was presented by NagaShima and Whitehead [Ref. 15] 
using dipole distributions. All three solutions were in 
close agreement. Sisto and Ni ({Ref. 20] developed a ‘time 
Marching" technigue applicable to an infinite cascade. Their 
results appear to be in good agreement with the finite 


cascade results. 


An analytic development, valid for low frequency blade 
motion, has been developed by Kurosaka [Ref. 11] and is 
currently being extended to the higher frequency range 
[Ref. 12]. These results are applied to an infinite cascade. 
Very recently, another analytic treatment has been presented 
by Verdon {Ref. 25] (see also Verdon and McCune [Ref. 24]}) 
in which an infinite series representation is used to 
express the influence of neighboring blades and wakes on the 
reference pasSage. However, serieS convergence was not 
always obtained, and some doubt exists as to its unifora 
validity Cher. 124]. A new analysis currently being 
developed by Verdon, applicable to the infinite cascade, 
does appear to yield uniformly good results [Ref. 26]. 
Actual wind tunnel tests of cascade flutter have been made 
by Snyder and Commerford [{ Ref. 19] and by Fleeter j{Ref. 7]. 


In the current study, a linearized method of 
characteristics procedure, based on a developement by Teipel 
inet. 213 for unsteady aerodynamics of « flat plate 


airfoil, is used to examine the pressure diStributions and 


ny 





4 


flutter boundaries for cascade blades having subsonic 
leading edge locus. It is an extension of previous work done 
by Chalkley (Ref. 5] and Brix [Ref. 2]. 


In addition, the method of characteristics was applied 
to the problem of supersonic flow through oscillating 
cylindrical shells to determine the pressure distribution 
and associated generalized forces acting on the cylinder 


walls. 


Although it is unreasonable to suppose that the 
two-dimensional cascade investigation can accurately model 
the complex mixed transonic-supersonic flow of an actual 
present day turbomachine, it is hoped that this analysis 
will provide a better understanding of the problems 
involved, and serve as a starting point for the flutter 
analysis during the design stage of modern high speed fans 


and compressors. 
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II. THEORY 


A. PROBLEM FORMULATION 


Consider a cascade composed of a semi-infinite array of 
thin, two-dimensional flat plate airfoils immersed in the 
supersonic flow of an inviscid, non-conducting, ideal gas 
having constant specific heats. Additionally, the flow 


field is assumed to be irrotational and isentropic. 





Figure 1. Cascade geometry 


The only restriction on the cascade is that the locus of 


blade leading edges must be subsonic; i.e., 


a (II- 1) 
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where @ is the stagger angle and a is the Mach angle 
1 

( sint!1— ). As shown in figure 1, the cascade is aligned in 
M 


the X~Y coordinate system so that the leading-edge of the 
first blade is at the origin and the blade itself lies along 


the X-axis. 


Perturbations in the flow field are caused by blade 
motion in pitch and plunge, assumed to be of small amplitude 
with simple harmonic motion. For a cascade with subsonic 
leading-edge locus, the blades are swept aft of the initial 
Mach wave. Thus, disturbances are allowed to propagate 
throughout the domain of influence (the region downstream of 
the Mach cone emanating from the leading edge of the first 
blade), and are not restricted to the region between blade 


passages. 
B. GOVERNING EQUATIONS 


Subject to the initial assumptions stated above, the 
governing eguations for this boundary value problem are the 


continuity eguation 


Se) = 0 II- 2 
> — +—— = = 
DT Ox ay’ ( 


Du 1g@p 

aout ae = 0 (21-32) 
DT pdx 

Dv 10 

et 0 (i= 3b) 
Dr poy 


the condition of irrotationality 
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—_ - — = Q (II- 4) 


er (i= 5) 


D ; ; . 
a denotes the material, or substantive, derivative of a 


fluid particle which is a function of position and time. It 


is composed of the local derivative (rate of change with 
respect to time) and the convective derivative (rate of 
change with respect to position) and is written as 

D( ) a( ) a 


a = ne + (Ve rC) (11=5a) 





where V is the velocity vector. The local sonic velocity is 
given by, 


d 
ce = 22) = bods (ek 6) 
Op s dp 
and 
P 
aie = constant (iE) 7) 


Taking the total differential of Eq. (7) 


ney 1 (y+1) 

ma) Cle — dp = 0 c= 25 

G) Pp yP@) Pp (il ) 
Rearranging 

dp p ; 

ae fee (II- 9) 

dp p 


Or 


PAI 





(II~10) 


C. BOUNDARY CONDITIONS 


The flow boundary condition (flow tangency) for a blade 
having an equation of the form 


Peet) =" 0 (II-11) 


is developed as follows. At each point of a solid-fluid 
surface, at every instant, the norawal component of the 
relative velocity between the solid and the fluid nust 
vanish. In cther words, the total time rate of change in 
F(X,Y,T) following a surface particie of the solid is zero. 


Mathematically, this statement is expressed as 


DF 
= = 0 (Ii~-12) 
DT 

For the cascade, the equations for the upper and lower 


surfaces of a blade are written 


Em(Xe ees) = Y = ¥ (X,T) = 0 (II-13) 
u u 
and 
ROG Y,1) =9¥ = Y (X,T) = 0 (II-14) 
L un 
Applying Eq. (12) to Eqs. (13) and (14) 
DF, OF OF, oF 
— = et it v= OO (II-15) 
DT OT ox oY 
or 
oY oY 
oe = 7 = 0 (aE 16) 
oT ox 
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. OF, 
Since —= 1. Likewise, 
Oy 
DF, Or eee 
DT oT ox 


(II~17) 


The other boundary condition that must be enforced is the 


Sommerfeld tadiation condition, i.e., 


propagate away from their sources. 


D. LINEARIZATION 


The normal velocity perturbation is written 


ead k, Linea 
’ ’ =a eae 
et ae )) at ox 
OY aY 

myer) = + 7 
get iS )) oT "ax 


From the assumption of thin blades and 


disturbances must 


(II-18) 


(II-19) 


amplitude 


oscillations, all flow quantities may be considered to be 


Small quantities which are linearly superimposed 


freestream guantities. Thus 


i= uu + u 


al 
u 
OQ 
+ 
Q 


Viv) 


while the density and pressure perturbations are 


Zp =p = a 
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onto the 


(II-20) 


(15-21) 


(II-22) 


(II-23) 





Ap = p-p (11-24) 


Also, because of the assumption of thin blades, Y and Y 
u 
are small quantities with respect to axial distances along 


the blade. Thus, substituting Eqs. (20) and (22) into Eqs. 
(18) and (19) and neglecting the effect of higher order 


terms, the normal velocity perturbation eguation becomes 


X,Y (X,T nf O%y II-25 
- - =— + ee —_ 
V(XeY eT) = 57 a (II-25) 
X,Y (X,f a ont eh OE 
Vv hi a on es a = 
CNG a = tN (II-26) 


To transfer the boundary condition to the axis, Eqs. 
(25) and (26) are expanded in a Taylor series about Y= 0. 


dv (X,Y=0*,T) 
nS 7 


v(X,Y (X,T)) = v(X,Y=0+,7) + ¥ 
u u oY 





1 @g2v(X,Y=0+t,7) 
Ss as ee ee + e@ © © (II-27) 
2 (bl (0) 2 


Ovi( 2-07, 1) 
ot (oy = v(X,Y=0-,T) + Yo —————————— + 


L oY 
1 a2@v (X,Y=0-,T) 
—Y2 —_________—_ + e ee (II-28) 
21. (OY) 2 
Neglecting higher order terns, 
wack, Yon T)) = vk, Y=0+,7T) (II-29) 
u 
V(XsY (eT) ) = Ul MeO 10) (II-30) 


and thus 
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i Ur oO = u JI-31 

v( oe ee ( ) 
oy, ay, 

Omg ky) = oscae aaa II-32 

ee) at a oax ( 


1 )9h 
For cascaded blades, it is assumed that the na blade 
st 
leads the motion of the (n-1) blade by an amount 6, the 


interblade phase angle. Expressing this and the assumption 
of sinusoidal motion, the equation for the upper and lower 
blade surfaces becomes 


ifwt + -1)8 
Y (X) te Teed (11-33) 
u 


Y (X,T) 
u 


ifw 4 
rmer” eae s (II-34) 


Yo (X,7 
eet) 


which reduces the flow tangency condition, Eqs. (31) and 
(a2), to 


OY if{(n-1)8} io 
e 


+ = LwY ——= - 

v (X, 0+,T) [ lw we + ee (11-35) 
oie aA: OY if{(n-1)8} if bs 

= = + — es 

v (X, 0-,T) [ lw et ) a” e (11-36) 
= iwT 
Now let Eg s a! = Y¥(X)e : Then the flow 
tangency condition for pitch is, 

Y (X) =- a (x = xo) (22-37) 


pitch 


where, 
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a= aie (11-38) 


while for the plunge, 


Y(X " =-heée II-393 
( Reiinide 5 (i ) 


Tar 
h = hie (II-40) 


Deflections are defined as positive downward, and angle of 
twist positive leading edge up. hy and a are dimensionless 
complex amplitudes, with xo the twist axis location (elastic 


axis) in percent chord ¢. Figure 2 shows the positive 
direction of all quantities used in defining the flow 


tangency condition. 


y 


undeflected 





Xo 


Figure 2. Positive @irections of Lift 
and Moment resulting from deflections 
in pitch (a) and/or plunge (h). 


To non-dimensionalize the flow tangency equations, 
distances are scaled to the blade chord and time to blade 
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chord divided by freestream velocity woe This leads to the 
introduction of the reduced frequency parameter 


WC 
k =— (II-41) 
u 


° 
which is twice the reduced frequency as defined by Garrick 
and Rubinow [Ref. 9] who scaled all quantities to blade 


semi-chord. 


Non-dimensionalizing Eqs. (35) and (36) and using the 
definition of reduced frequency above, the flow tangency 


conditions become ° 
pitch: 
Wilk, Ot, c) ikt 
—————— = - a {cos[ (n~1)8] - k(x-x_)sin[{ (n-1) 8 J}e 
uy 0 ce) 
ikt 
_ da) {ksin{ (n-1) 8] + k (x-x } cos[ (n-1) 8 Be 
(II-42) 
plunge: 
Valk, 0 t,t) ikt 
——— = h) {ksin{ (n-1) 8] - ikcos[ (n-1) §J}e 
°o 
(II-43) 


To linearize the expression for the pressure, Eqs. (10), 
(20) and (24) are used. Expand og. (10) in terms of 
perturbation quantities. Through Eq. (48), second order 


terms are neglected whenever they occur. 








(p+ Ap) 
c2 + 2cc = = — II-44 
2 0 (1 +42) ( } 
p 
° 
Ap , : ; 
As p <1, the denominator can be expanded in a geometric 
is] 
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series of the forn 





1 
=41-et e@- eX tooe (II-45) 
1 +e 
to yield, 
Ap 
eemt 2CG = y(o. + 4p) (1 = —-) (II-46) 
to) 0 (3) i. 
Expanding, 
c2A 
c? + acc. = (1 Bye Np (II-47) 


y Op 


However, as Ap and Ap are both small quantities 


(II-48) 


> 
ro 
Pu 
© 
Oo! 
ft 


Expanding Eq. (47) , substituting into Eq. (46) and 


rearranging terns, 





Ap 2 2ye =E 
pe y-1 0 (y~1) ¢, 
2yc : , 

—— CC 7 Using Eg. (43) to expand the term in 

(y~ 1) ¢, 
brackets yields the linearized form of the pressure 
difference 

: ro 50 
0 Dao ae ae (II-50) 


' 28 





TIL. ) SeLuOD OF (CHARACTERISTICS 


A. CHARACTERISTIC DIRECTIONS 


To introduce the method of characteristics, it is 
important to recall that for supersonic flow there are 
distinct regions of influence and dependence. Hence, curves 
must exist along which the vaiues of the flow variables and 
their derivatives are known, yet the flow field in the 
immediate neighborhood of these curves is indeterminate. 
Therefore, to obtain the directions in the x-y plane of 
these curves, which are called characteristic directions, oL 
more simply, characteristics, one must prescribe that the 


derivatives across these curves are indeterminate. 


Mathematically, this means that the dependent variables 
must be continuous, but their first, or higher order, 
derivatives may be discontinuous. Physically, the 
characteristics act as information carriers, in that ail 
disturbances propagate along then. In particular, they 
represent the paths of Mach waves, and are therefore 


referred to as Mach lines. 


fo facilitate the analysis, it is desirable to introduce 
new coordinates € and 7». fhe governing equations are 
written in terms of the two arbitrary, but intersecting 


curves 


my 
{ 


ee kr, y) constant (III~1a) 


and 


n (x,y) constant (III-1b) 


3 
it 


However, the equations are more Suitably transformed to the 


arbitrary coordinates if the dependent variables are UW, v 
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and Cc. Here, air is assumed to he a perfect gas with 


constant specific heats c (constant pressure) and c 
* Pp Vv 


(constant volume). Rewriting Eq. (II-10) as 


G2 = yRO (III- 2) 
Ro= ¢ = ¢ (III- 3) 
p v ; 
where R is the gas constant for air, and 8 is 
temperature 7 degrees Rankine. Taking the total 


differential of Eq. (2) and dividing by C2, 


dc as 
= = — (TII- 4) 
c 8 


Taking the total differential of Eg. (II-10), 
1 dp 


2éde = y-dp - p— III~- 5 
aan) Pa: ( ) 


oc, dividing by C2, 
2— =—- — (III- 6) 


Now from the First Law of Thermodynamics and the fact that 


the flow is homentropic, 


P 
@ds = du. - —dp (III- 7) 
int pe 
Horeean ideal gas, du = c dO and Eq. (7) may be written 
int Vv 
as 
dp 
eds = c d@ - RO— (III- 8) 
yf P 
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Dividing by c 9 
v 


1 de dp 
is 0a = (y= i) (1ii- 9) 
Cy 8 Pp 


and using Eqs. (4) and (6) yields 


1 dé 1 
—ds = 2— - (y-1)—dp (III-10) 
é P 


ey 


But because the thermodynamic properties of a continuous 
flow field are to be considered as continuous, single-valued 


scalar point functions, Eq. (10) must hold for any change of 


the fluid properties s, ¢ and p. Hence, it may be written in 
terms of material derivatives. Then, subject to Eq. (il—3)e, 


Eq. (10) becomes 


1 Ds De 1Dp 
—— = 2—-— - (y-1)-—— = 0 (ETi-—11) 
c Dr CDT pbtT 


Then, uSing Eq. (Ii-10), Eq. (11) reduces to 


— =—— t— (III-12) 


Now the perturbation equations are introduced. Writing Eq. 
(II-9) as 


Dp _ Dp 
— - o— (2IT=13) 
DT DT 


and substituting into Eg. (12) allows Eq. (II-2) to be 
expressed as 
2 ac 2 ac ~ 6u dv 


Stl ee ee TIrI-14 
y-1dT y-1 AX 0dXx OoY 


are 


heglecting the effect of higher order terms. This is the 
form of the continuity equation which is to be used in the 


transformation. 
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Moecencerethe first Euler equation, Eq. (11-3), in a 
form suitable for transformation, use Eq. (11) in 


differential forn 


an 1 
at = (y-1)—dp (II1I-15) 
c P 


With Eq. (II-9) this yields 
ede (I1I-16) 


Expanding the total derivatives and equating the partial 


derivatives, one obtains 


——¢-— {Tii-17) 


ee ys = = (III-18) 


Introducing the small perturbation quantities and neglecting 
second order quantities yields the desired form of the first 
Euler equation, 

ou ou 2 dc 


— + —— —— = Q III~19 
eae ” oax y-1° 00x ( 


— 


== = as (III-20) 


is used in place of the second Euler eguation as one of the 
governing equations to be transformed. Then the continuity 
equation, the first Euler equation and the condition of 


irrotationality form a system of three equations in the 


three unknowns U, v and cC¢. For the assumption of simple 


32 





harmonic motion, the Teipel amplitude functions [Ref. 21] 


are introduced, 


iwT u- iu 








U(X,Y)e = ° (II I-21a) 
u 
Q 
V (X,Y a wee III-21b 
e = — ~- 
(X,¥) 7 a, 
iwT A | es cy 
CXx Ve = —~(—-)——___$ (ILI-21c) 


y-1 He Cc, 


where U, V and C are complex non-dimensional amplitudes. U 
and V are the perturbation velocities and C is the 
perturbation in sonic velocity which will be related to the 


pressure coefficient through Eq. (II-49). 


Suostituting Eqs. (21) into Eqs. (14}, (19) and (20), 


and then non-dimensSionalizing as before yields, 


du Ov ae 
— + /m2-7 — + M2— + ikm2@c = 0 (IL I-22a) 
ox oy Ox 
ou oc : 
— + —+ ikuU = 0 (II I-22b) 
Ox Ox 
ou ov 
—- /M2-1 —= (III-22c) 
oy x 
For € = constant or 7 = constant, both ¢ (x,y(x)) and 
m {X,Y (X))} are implicit functions of x. Hence, 
dg = ae +€ dy =0 ffl 1-234) 
Y 


and 
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dy = 7 dx + m dy = 0 (ILTI-23b) 
x Bg 
Therefore, 


sts - — (III- 24) 
y ‘ly 
and the labeling of the characteristic directions is 
arbitrary. By convention the left-running Mach line is 
normally denoted € and the right-running Mach line, n° 


To solve for the characteristic directions in the x-y 
plane, Eqs. (22) are rewritten along lines of constant € 


and 7. 
eaUeeeteyte- teen Vo etence C ==> U - Me, C = 
x é yé x € Vy n 1s n 
M2-1 y V — ikM2c (III-25a) 
Y 7 
Ut C =f Y= (shill III-25b 
ae S é 7x Ty n 
€u - Mi2-1 € v =- 7 U + /h2-1 » V (III-25c) 
y. x ey xn 


Now Cramer's rule is applied to solve for the 
derivatives of the dependent variables with respect to €& 


along 7 = constant. To satisfy the condition that these 


derivatives be finite requires that each be of the fora ie 

es » indeterminate. Equating the determinant of the 
e « e c ; 

coefficient matrix of {U,C,V } to zero yields, 


cele (== 1) = <2 ] = 0 (III- 26) 
x x y 
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to which there are three real solutions, each of which 


defines a characteristic direction in the xry plane. These 





are, 
dy : 
—=-+f =0 (III-27a) 
dx x . 
ay on 1 
mae fk III-27b) 
ax ey ¥M2-1 ! 


The same results are obtained when one solves for the first 


derivatives with respect to 7 along lines of € = constant. 


EG (27a) gives the streamline slope in the flow field, 
while Eq. (27b) gives the slope of the two curves € = 
constant and 7 += constant. It should be noted that these 
are the same characteristic directions as obtained for 


steady, two-dimensional flow. By convention, 








dy 1 
ae = + IIli-284 
ae =const /M2—-1 ! 
dy 1 
(e— = = (IITI~28b) 
dx m=const Me-1 
dy 
— Se of (III-28c) 
dx strean 
For an arbitrary function f£ = f(x,y), the total 
derivative is 
are of 
df = —dx + —dy (LILI 29) 
ex Oy 


For £(x,y) taken along a curve as defined by Eqs. (1), 
however, y is an implicit function of x and Eq. (29) is 


rewritten as 
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df _ df | dfay 








—=s— (III- 30) 
ax bse dyax 
Therefore, along a streamline, 
df Of 
— = (III- 31) 
dx stream 0x 
while along the Mach lines, 
df Of 1 of 
en) = oat ee: (II I-32a) 
dx €=const Ox Jfi2-1 dy 
af _ of Toot (Ir 1-32b) 
ax n=const Jim /N==1 oy 


Then, solving for the partial derivatives in terms of the 


ordinary derivatives, 


Of df ae 
—_ = — ItIi- 

Ox Uae ate ( a 
Be = c + ue IIL Ie 3}.33}8) 
ax st = =const tp ccenait , ( 
he = = ie = ILM 3S 
oy 7 TUS, =const econ ( - 


The compatibility relations are obtained by writing Eqs. 
(22) using Eqs. (33). Thus, 


i du du ae av 
= =o + eee rs — - — 
5 (Re le J (He -1)\[ os ~ aes J 
H2 ac 
ade tale te LMSC = 0 (II I-34a) 
xe dx 7 
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1 dU du ; 1. dc 
ae —— & f= anff ff—e 
at f ) eee pie) 


dice E 
(IT I-34b) 
= a - [ a + a iG= 0 (II I-34c) 
eee Gale * G2), 


Eq. (34b) is multiplied by M2 and then subtracted from Eq. 
(34a). The resulting equation is first added to, and then 
subtracted from Eg. (34c) to give two of the compatibility 








equations 
as a ik = uU-C 0 III-35a) 
= - (— +a = = =355 
le tre oon ) ( 
au dv _ M2 
gant) eto Ges) YR (Use) = 0 (II I-35b) 
ax 7 ax n 2 


The third compatibility relation is obtained from writing 
Eq. (34b) along a streamline. 
du dc 


(-—) (=) + ikuU 
dx stream dx stream 


" 
io) 


(II I~35c) 


The governing equations have now been reduced to a system of 
three ordinary differential equations which are implicitly a 
function of x, with U , V and C complex amplitudes. fhus, 
Eqs. (35) are really six equations, with three real and 
three imaginary components. 


To write the boundary conditions in terms of the Teipel 
amplitude functions, note that 


ikt 
v(x,0,t) = vu Mi2-1 V(x,y)e (III- 36) 


from the definition of the amplitude function, and therefore 


the tlow tangency condition is, 


a7 





ip v(x, 0,t) 





V= V(x = Ei i=_37 
( rY) Juen~t uy, ( ) 
v(x,0,t) . ; ; 
where ——————_ is given by (11-40) for pitch and (II-41) for 
u 
° 
plunge. 


B. FINITE DIFFERENCE EQUATIONS 


To write Eqs. (II1I-34) in finite difference form, the 


computational nolecule shown in figure 3 is used. 








current Mach 
line 


previous Mach 
line 





Figure 3. Computational molecule for a 
general fiow field point. 


If F is an arbitrary complex flow quantity, derivatives 
with respect to x along the characteristics, in finite 


difference fcrw, are written, 


dF Faas. F 
(=) = -<20 |! ., (III-38a) 
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2a es l2 (III-38b) 
“az é Lax 
aF F a IP 
(—) = oe 2 (IL I-38c) 
dx n 5 AX 


The values of the perturbation guantities in Eqs. (38) 


are averages, so that 


F =-—(F ap fe HID Be 
‘ ae 2! 11 22 ‘ ay 
(F) 2, F + F ) (ITI+39b) 
Eo 42 22 
F : F + F III-39 
= = ad (e 
ee 5, ( ) 


ieeceneral Flow Field Equations 


With these relations substituted into Eqs. (35), the systen 


of equations becomes, 


U + C - AU ouK (III-40a) 
22R 22R i 220 12R 

U + C + AU = K (II I-40b) 
221 221 I 22R 121 

U - Vv -~B(U - Cc ) ae (III-40c) 
22R 22R I. 228 221 34R 

U - Vv +B (U - Cc )=K (ILI-404d) 
221 22 1) 22R 22R 341 

U + V - B(U -C )=K (III-40e) 
22R 22R D221 220 56R 
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U + V + Bec U -c ya is (II I-40£) 

221 221 I 22R 22R 561 
where 

K =) + C + AU (IZI-40g) 
12R 11R 11R (ae 

K = + Cc - AU (II I-46h) 
121 111 111 I 11R 

K = 9 - Vv +B (U Pere.) (II I-40i) 
34R 12R 12R eet 27: 121 

K = 9 - Vv -B(U - Cc ) (III-404) 
341 12 121 Tee OR 12R 

K = + Vv + B (U - Cc ) (III-40k) 
56R 21R 21R Tee 20r 211 

K =) + V - B(U -Cc) (II I-401) 
561 ase 21 feet R 21R 

and 
1 

A = —-kAx (III-4 1a) 

2 
1 Ie : 
B = ~k(——j Ax (TII-41b) 


I 4 M2-1 


As shown by Teipel [Ref. 21], solving these equations gives 


: al 1 
ee A,8,)[$( Kz an*Ksea)- 81X19) | +28, H K5.4)*Kag))*8; Kin 


22R _ 2 2 
1-A.B 
( 8, + (28)) (LII- 42a) 
7 iF [3 | 
; i re outenre Kee 6 Kio. | - 26d o(K. 06K. )-BrKi 
221 (1-A,B))? + (28,7 
(1II-42b) 
1 
= -—( K - K ) (IIi-42c) 
22R 2  S6R 34R 
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These are 


located in 


For 


K A U 
121 221 E225 


(IX 1-424) 


(II1-42e) 


(II I-42£) 


the finite difference equations for a grid point 


the general flow field. 


a grid point on top 


of a blade, the 


computational molecule of figure 3 is modified as shown in 


figure 4. 


Figure 4. 





Computational molecule for a grid 


point on the upper surface of a blade. 


Perturbation quantities at grid point a are considered 


zero and the normal velocity ue 


Sant F 
Pp ag 


(pitch or plunge at grid 


is prescribed by the blade movement as given by 
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Eq. 


as given by the flow tangency condition. Blade 


(11-37). 


U c + 
221 221 


U B U + 
22R i BAN 


+ 3 


U U 
221 ee 2k 


K and K 
12 56 


K 
34T 221 


BC 
dee 2 


C 
I 22R 
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are as before, 


Here the applicable eyuations are, 


(IT I-43a) 
(IL1I-43b) 
(II I~43c) 
34R 
K (ILI-434d) 
341 
but 
(ILI~43e) 


(II I-43£) 


position is 


always given relative to the leading edge of each blade. 
Solving, 
1-A 8) [k eye | : 
= Pe oon misaReT ey) il sel 341" pier 
22R (1-A,B))? + (2B, )? 
(111-444) 
1-4, “Kear 8a: “ek 
fee Pr) LXser —Ksay* yin 28, K34R BK] 
22] (1-A,B,)? + (2B,)? 
(III -44b) 
Cc = 7K - + U (ITI-44c) 
22R 12R 22R q 221 
c = K - - AU (III-444) 
221 Ueade 221 I 22R 
and V are known. 
22 





At a grid point on the lower surface of a blade, the 


normal perturbation velocity at grid point ES is also 


known. The computational molecule used is shown in figure 


Di 





xX 


Figure 5. Computational molecule used for a 
grid point on the lower surface of a blade. 


For this grid point, the finite difference eguations 


become, 
U + C - AU =k (II I-45a) 
22k 228 E221 a2k 
U + C + AU ah (LII-45b) 
221 228 i 228 127 
U ~ B(U - Cc a ak + K (III-45c) 
22k ig 220 Do 56R 34R 
+B ( U - Cc i) eS Is + K (IT 1I-454) 
221 I Z2R 22k 561 341 


where M3 is the same as Eqs. (39a) and (39b), but 
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K = U -Vt ¢ 3G -c ) (ILI~46a) 


K = U - Vv - B (U =" Clw| (ILI-46b) 
561 121 121 I 12R 12R 


K = V (III-46c) 
34R 22R 


K = V (ILI-464) 


Solving for the perturbation quantities at grid point 22 


yields, 
4 _1-4,8)) [Keeat S3an7 21h ] +28 keg Pee | 
aa (1-A 8, )%+ (28,)? 
(111-47a) 
7. _1-8,8) [Kyor + Kar * 8rX ron) 281lsent K3an 8rKro| 
pel (1- A,B, )% + (28,) 
(IL1-47b) 
Gel aK =u +AU (ITI-47c) 
Boren 2h = 22B =~ 221 
Oe= Ki =U - AU (III-474) 
eo 127) 220° t 22k 


and V are known. 
22 


4. Initial Left Running Mach Line 


To obtain the equations for the initial left running 
Mach line (€ = constant), the computational nolecule of 
figure 5 is used. However, it is assumed that grid points 
is and on are just in the freestream and the limiting 


condition x--0 is taken. For x-+0, AS and ae Thus Eqs. 


ay 





(38) reduce to, 


i 1+ Vv =O (II I-48a) 
22 20 
Jeet c = 0 (II I-48b) 
DD Da 
OL, 
= ="y == c (III- 49) 


Now Eq. (35a) becomes 
M 
(—). = - ik—J {EL 1L= 50) 
x M 


which can be integrated along € = constant. Or, 


Me 





U = 2 exp[-ik x] (L1iI=-"571) 
x=0 


42-1 


U F is known from Eqs. (49) and (37) Thus, 
x= 


2 M2 
x) - V 0) sin(k 
) oon ) ( 


a 
it) 








av (0) cos (k x) 
R 


22R 22 M2-1 H2—t 


(III-52a) 


U 








V (0) sin (k 
R 


M2 
U x) = V (0) cos (k 
221 22 2 dé 


M2 
x) 
jy cee | ae M2-1 


(II I-52b) 


with the appropriate form of Eq. (37) used for a Then, 


V =-U (II I-52c) 
22R 22R 
=-U (III-524) 
Dog 220 
Cc a (II I~-52e) 
22R 22R 
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=- U (III-52£) 


Eqs. (52) are also used to compute the perturbation 
quantities at the leading edge of the upper surface of the 
first blade. 


5. Initial Right Running Mach Line 


The developement for the perturbation quantities 
along the initial right running Mach line (= constant) is 


analagous to the development for the left. Here, 26 and Poe 


of figure 2 are assumed to be just in the freestream. Eqs. 
(38) reduce to, 


i ae) (III-53a) 
20 22 
U.. ee, = 0 (II I-53b) 
22 22 
OL, 
U=Ve=-C (III- 54) 


Then Eg. (35b) becomes, 





x (III- 55) 


Integrating along n = constant yields, 











M2 
U = U exp[{-ik x III- 56 
o Pl M21 | ( 
OL, 
M2 M2 
U = V 0) cos (k x) + V 0) sin (k x 
22R 20R | ) ( M2-1 ve ( M2-1 


(ET 1-5/4) 
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a 


=- V 0)sin(k 
22R° ( 


M2 Me 
x) + V 0) cos(k 
M2-1 221° ) ( 








x) 


U 
22u M2-1 


(IL1I-57b) 
Again, the appropriate form of Eq. (37) is used to define 


VV (0). Then, 


22 

v = U (III-57c) 
22R 22R 

V = JU (ILI-574) 
per 221 

C =- U (ILI-57e) 
22R 22R 

C are (ILI-57f£) 
22h 221 


Eqs. (57) are used to compute the perturbation guantities at 
the first grid point on the lower surface of the first 
blade. In the computer program, lower surface points at 


grid point a are subscripted 33 or 55 for pitch and plunge 


respectively. 


6. First Grid Point on a New Blade 


For the first grid point encountered on any blade 

other than the first, the computational molecule of figure 6 

is used, The new blade is assumed to protrude an 

infinitisimal distance past grid point ee which is labeled 

FE and F in the figure. Then, the equations developed 
22U 2k 


for grid points on the upper and lower surface of a blade 


are used to compute the perturbation quantities at aon and 
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F respectively. 
22L : y 


y 





Figure 6. Computational molecule used for 
the first grid point on a new blade. 





Figure 7. Computational molecule used for 
a grid point in the wake. 


It is assumed that the first wake grid point is an 


Infinitisinal distance aft of the blade. To calculate the 
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flow field quantities in the wake, the computational 
molecule is divided into two triangular molecules as_ shown. 
Then, the equations developed for a grid point on the upper 
surface of a blade apply to Ee while those for the lower 


surface of a blade apply to F 5 The computational 


221 
molecule is shown in figure 7. 
U #c - AU = (III-58a) 
22RU 22RU L 220 12RU 
U + C +AU = ik (III-58b) 
2210 221U I 22RU 121IU 
U + Vv +B (C - U hesok (III-58c) 
22RU 22RU 1g 221U 221U 56R 
U + V - B (C - U j=) kK (III-58d) 
2210 221U I 22RU 22RU 561 
and, 
U o - AU = K (IITI-59a) 
22RL 22RL I 221IL 12RL ; : 
U a & +AU an (II I~59b) 
220 22h ie 22 RL 12IL 
U - Vv + B(C - U j= ok (III-59c) 
Z2nL 22RL i 221IL 2204 34R 
U - Vv = 8B (€ - U \ SR (111-594) 
220 221k I 22RL 22RL 341 
where, 
K = ii + C +AU (II I-60a) 
12RU 1TRU 11TRU I i1LU 
K = JU Ae -aJvU (IL1-60b) 
T21U 11LU 1110 I 11RU 


49 





K = U iG + AU (II1I-60c) 
12RL 11RL V1RL © Piih 

K = § - AU (IZI-604) 
1270 MMIL TEL I 11RL 

K = U - Vv + B ( U - Cc ) (II 1I-60e) 
34R 12K 12R I 121 122 

K = U - Vv - B(U - C ) (ILI-60£) 
341 ea 121 ug 12R 12R 

K = U + V + B ( U - Cc ) (II 1-609) 
56R 21R 21R i 211 211 

K = 9 + V - B(U - Cc ) (III76Gh) 
561 211 2A E I 21R 21R 


Egs. (58) - (60) form a system of eight equations 
with twelve unknowns. However, in the linearized problen, 
the wake cannot support a pressure difference, and the 
normal velocity across the wake must be a constant. AS a 
result of these two conditions, to first order effect, there 
will be no discontinuities reflected from the wake. 


Therefore, the relations, 


V =v =v (III-61a) 
22RU 22RL 22R 

v = Vv = Vv (III~61b) 
221U 22IL 221 

é =a Sc (ILI~61c) 
22RU 22RL 22R 


Cc Cc =C (ILI=614) 
2270 221L 221 


are used to reduce Eqs (58) and (59) to a determinant eight 


by eight system of equations which are: 
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a 


U + C - AU = K (ILI-62a) 
22RU 22R I 2210 12RU 

U + C +AU = K (III-62b) 
2210 221 I 22RU 1210 

U + V a BAG = {i ) = K (III~-62c) 
22RU 22R I 221 2210 56R 

U + V Dae ee a ) = K (ILI-62a) 
22RU 221 I 22R 22RU 561 

and 

U +c - AU = K (IIT I-63a) 
22RL 22R I 221L 12RL 

U 2 + AU = K (IL I-63b) 
22IL 221 I 22RL 12IL 

U + V wedee (IC - U ) = K (ILI-63c) 
22RL 22R I 221 22IL 34R 

U + V he C - U ) = Kk (II I-63d) 
22IL 221 I 22R 22RL 341 


Egs. (60) remain unchanged. 
C. THE PRESSURE COEFFICIENT 


The non-dimensional pressure distribution along the 
surface of a blade is determined using the Linearized 
pressure difference of Eq. (I1-48) and Teipel's sonic 
velocity perturbation amplitude function, Ege (21c). 
Rewriting Eq. (II-48) as, 


2 e 
eee pce (-———_—_——_} (III- 64) 


and then dividing by the freestream dynamic pressure yields 
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® 


Daeee PRE 
Les = 2C(x,y)e (1ii- 65} 
2 


Hence, the ccmplex pressure difference across a bliade at 


grid point se resulting from oscillation in pitch or plunge 
is, 


Pp -p ikt 
POM oe) = Saat Se Te pe (II I-66a) 
3YpM2 22L 22U 


This expression for the pressure difference across a 
blade appears to be independent of interblade phase angle 6. 
However, the interblade phase angle was used for the 
computation of the normal velocity through the flow tangency 


equations, and they were used in the determination of Ce 
Therefore, to compute the pressure distribution along the 


th 
n blade when it is in a specific position reguires that 
the position be specified in terms of the position of the 
first blade. In general, the complex pressure is 
; ikt 
P(x,y,t) = oes + ae (III- 67) 


th 
with and =e defined by Eq. (66). The position of the an 


blade is 
(kt) = (kt) + (n-1)8 (III- 68) 
n 


Then expanding Eq. (67), 
P(x,y,t) = [P cos(kt) - P sin(kt) J 
R I 


PUP E Sot tke) + eeceet kt) | (III- 69) 
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Figure 8 shows a complete oscillation cycle for a_ blade 
undergoing simple hamonic motion in pitch. The extension to 
plunge is cbvious, and the following discussion is 


applicable there also. 





Figure 8. Blade oscillation cycle. 


The two conditions of primary interest here are when the 
h 
n blade is in the initial position 


(kt) = 0 (III- 70) 
n 


or in the mean position, pitching leading edge up, 


(kt) = 270 (IZI- 71) 
n 


Then solving for the position of the first blade when the 
n : blade is in the initial position gives, 

(kt) = - (n~1)8 (III- 72) 
and in the mean position, 


(kt) = 270 - (n-1)8 (III- 73) 


Except for flutter analysis, the real component of the 
pressure is normally of primary interest, and can be found 


by solving Fq. (69) with (kt) expressed as a function of 
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(kt) . For ease of showing the relation of the real 
n 


component of pressure for the two blade positions above, 


consider the special case of 8 = 0. Then, 
Re {P(x t = Pp Dib ia 
(P(xeyrt)} = PS (III-74a) 
Re {P(x t = Pp ILII-74p 
{P(XsY, ac I ( ) 


and thus, 


Re{P(x,y,t = Im{P(x,y,t ILI-75a 
(Px AF Some ee ( 


Re {P(x t = Im{P(x,y,t III-75b 
e(P(X,y, 3 2070 {P (X,Y, Jee ( ) 


Having obtained the pressure distribution over the 
blades, the dimensionless lift per unit span and the moment 


per unit span sguare are found fron, 


1 
L(x,y,t) = = feccy.tax (ILI-76a) 
0 
and 
{ 
M(x,y,t) = = [xg P ceore tax (L1I-76a) 
° 


where L(x,y,t) and M(x,y,t) are complex quantities which 
include the effects of blade oscillation in pitch and 
plunge. Lift is positive upward and moments are positive 


Clockwise (see figure 2). 
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IV. FLUTTER ANALYSIS 
A. EQUATIONS OF MOTION FOR A SINGLE BLADE 


Lane [Ref. 14] showed that cascade flutter could be 
analyzed by considering only a single oscillating blade. 
Following Fung {Ref. 8], the equations of motion for a flat 
plate airfoil immersed in an inviscid flow are derived from 
a force balance of the inertia, elastic and aerodynamic 
forces. For bending motion, h is meaSured positive down from 
the elastic axis. The torsional motion is poSitive clockwise 
about the elastic axis (refer to figure 2). It is assuned 


that a pair of springs, with spring constants Ko and ore 


oppose the blade motion. 


The total inertia force per unit Span is, 


F =~ (4 h + S %) (IV- 1) 
i b b 


and this force exerts a moment about the elastic axis s eye 
M =- (I @ + Sh) (IV- 2) 
a b 


where 

a ~ total blade mass per unit span 

oe - wing static moment (about x4) 

Io - wing mass moment of inertia (about x) 
The springs produce restoring forces and moments of 


F = - hK (IV-3a) 
ie h 
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Bie dhs (IV¥-3b) 
rt 


Then the equations of motion for the blade are, 


Lift/span = Mh + S @ + hK (IV-4a) 
b b h 


Moment/span = De + I 2 + ak (IV-4b) 


The assumption of simple harmonic motion is expressed as, 


iwT 
h = hie (I1V-5a) 
iat 
a= ae (1V-5b) 
where b and a are complex amplitudes. The blades 


oscillate with a constant interblade phase angle §&. A phase 


difference exists between twist and flexure also, but cannot 
be determined until a flutter solution is obtained. Then, 
the ratio of blade deflection amplitudes will give this 


phase difference. 


With Eqs. (5), the equations of motion for the blade 
become, 


5 = - 2 - 2 x -6 
Lift/span w oe w eecG + a (1V-6a) 


Moment/span = - w2S h -w?Ia +aK (IV~-6b) 
b 0 a 9 o a 


The natural frequencies of twist and flexure are found by 
assuming that each mode can exist separately, in the absence 
of exciting forces. Then, 


Ww = JM (EV — 74) 
b 
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cae he (IV-7b) 
a. a a 
OL, 
K = Mw2 (IV-8a) 
h h 
= 2 - 
Ky Se (IV~8b) 


and therefore, 


iw 
i = _— 2 ~” 2 =_ 2 — 
Lift/span [ how an al: aw ay the (IV-9a) 
iwt 
Moment/span = [-hw?2S - a Ww?2I - I w#?) je (IV-9b) 
9 b 0 a aa 


B. TWO DEGREE OF FREEDOM FLUTTER EQUATIONS 


The method used to investigate flutter in this report 


follows the procedure and notation used by Garrick and 


Rubinow [Ref. 9], who derived the lift force and moment 
equations, 
i 4 pb3w2 ees L nl z rab 
= e —2 +i + oe! 
G/R hie iso Fs a y) J 
(LV-10a) 
M 4 pb4w2 yet De cM + iM) (M+ iM )] 
G/R ‘aaa Sed al ae 4 
(IV-10b) 
where L was defined as positive down. Equating the 
G/R 


magnitudes of the lift from Eqs. (9a) and (10a), one 


obtains, 


oe] 





Bog iL eek) ot a (Lt LL ) 9) 
— to = a 1 = = 
Sue Qs i‘ 3 Bee 
(IV-11) 


and by equating the amplitude of the moments, 


h, 
—(M + iM - px +a (M + iM - 2¢90 xX) = 0 
a a ae 0M; i) aaa 








b 
(IV- 12) 
where the non-dimensional terms are defined as, 
Ss + are - ee ‘= real and imaginary component 
of lift and moment due to plunge. 
(ieee 1b), (sh) 6 6+ cil } - real and imaginary component 
3 4 3 ur 
of lift and moment due to pitch. 
: My 
pe. - wing density parameter ( ) 
4 pb2 
i - ratio of natural frequencies of flexure 
wh 2 
and torsion LG) 
Q 
Q - j 2. 
& the quantity ee 
x - location of the wing center of gravity measured 
a 
Sp 
from the elastic axis (——) 
Mpb 
: : : Ty 
t + radius of gyration of the wing ( ) 
a M,b2 
X - ratio of natural frequency in pitch to the 
Ww 2 
circular frequency () 
Experimental results indicate the Ww > w in modern 
a 


turbomachines [Ref. 18]. 
; h 
Eqs. (11) and (12) form an eigenvalue problem in = 
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and (a) which has a non-trivial solution if and only if the 


determinant of the coefficient matrix vanishes. FEquating 
the determinant to zero yields a complex function in the 
unknown X, each component of which must be identically 


equal to zero. These equations are, 


2 Q x2 Q Lb - +Q (M - Q)]X + C = 0 (LvV-13 
ee + [ me ; ) ah A a3 = ; ( a) 


(QL +4¢02M4 )X¥ +c = 0 (IV¥-13b) 
a2 h 4 I 
C= yw + L - (mM -2 -Lr2- 2) +D 
ea Be ee Gg? eee a R 
(1V-13c) 

on it + - - M -Lr2]) + D IV-13d 
I aa mas 2 u! 4 2 = I : : 
p= LM = L-§. = LM + LS (IV-13e) 
R 3 31 gan 4“ 2 

D=LM ~-LMN +LM -L4S (LV-13£) 
I 14 uf 2°3 3 2 


For a fixed cascade 
geometry and interblade 
phase angle, Eqs. (13) 
are all transcendental 
functions of the reduced 
frequency Ky OF (17k); 
more accurately. The 
solution of Eqs. (13a) 





and (13b) is found using 
the Theodorsen method. Figure 9. Solution plane for 


the Theodorsen method. 
/Re {X} and /1m {x} are 
tabulated and plotted versus (1/k). The flutter point is 


defined by the intersection of the real and iwaginary curves 
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Bi YX « The approximate shape of these curves is Shown in 
figure 9. Since (1/k) appears transcendentally in every 
term of Eqs. (13), the curve of /Re{x} is not exactly 
parabolic, and for /Im{xX}, not actually linear. 


At the flutter point, the non-dimensional flutter 


frequency is, 


we 1 
ay (I V- 14a) 
a 4X 
and the non-dimensional flutter speed is, 
Ur 11 
SR (IV-14b) 
we kyS 


k is defined by Eq. (II-39). 
C. SINGLE DEGKEE OF FREEDOM FLUTTER 


. For single degree of freedom oscillations in pitch, the 


equation of motion for the blade is, 


Moment/span = I a@ + aK (IV- 15) 
Qa. a 


which with the assumption of simple harmonic motion becomes, 


lwf 
Moment/span = [- w2I a_ + K q_Je (IV- 16) 
a 9 a 9 


The mechanical damping of the system is introduced by 
assuming the spring opposing the displacement in pitch to 


Maven a Spring constant of the form K (1 + ig) [Ref. 8]. 
a a 
The natural frequency in pitch is determined as_ before. 


Then, equating the amplitude of this expression for the 
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h 
moment with that of Eq. (10a) with “nd = 0, one obtains, 


Oty =9) +e = 0 (1V-17a) 
a 3 


M +g X = 0 (I V-17b) 
4 a a ; 


The solution method is the same as for the two degree of 
freedom problem, and the non-dimensional flutter frequency 
and speed are found from Eqs. (14). A similar development 
is required to obtain the equations for single degree of 
freedom plunge. However, based on the results of this, and 


other works, flunge oscillations are always stable. 


For both single degree of freedom and coupled 


flexure-torsion motion, the flutter calculations for the n 
blade must be repeated for varying interblade phase angles 
until the minimum flutter speed is found. This speed is 
termed the critical flutter speed. For all speeds below it, 
cascade operation will be stable. 


An alternate method for determining the onset of flutter 
is to compute the real component (the in phase component) of 
the work per cycle of the blade on the freestrean, 

2M 
Work/cycle = Re{/ (HM &dt + L hdat)} (IV- 18) 
0 . 2 

As long as the blade does positive work on the freestrean, 
the cascade will be stable. However, if the aerodynamic 
damping decreases to the point that the airstream does work 
on the blade, flutter will occur. For single degree of 
freedom pitch oscillations, with zero structural damping, 
the sign of the imaginary component of the moment determines 
whether the work/cycle is positive or ne ,.tive, and car 


therefore be used to indicate the onset of flutter 
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[Ref. 4,18). 
D. DEFINITICN OF TERMS 


The complex pressure difference at each grid point along 
a blade was shown by Eq. (I11-66) to be a direct function of 
the sonic velocity perturbation amplitude. Eqs. (III-76) 
then gave the non-dimensional lift force and moment per unit 
Span. These integrals were numerically computed using a 
trapezoidal integration routine. If F represents the 


integrated value of (m+1) discret values of £f(i), then, 


1 
P= f(x) dx = (5 (£Q) + £(m+1)] + > £ (i) Ax (IV- 19) 
i 


To relate the non-dimensional expressions for lift force 
and moment as derived by Garrick and Rubinow to _ those 
computed using the Teipel amplitude functions, the following 


quantities are defined: 


in). liwT 
P(@szy,t) = Cx Sa = one) Je (IV-20a) 


Then the dimensionsal lift force and moment per unit Span 
are, 


Piet = ; L( its = : = Ib L 
1 =puec x —puec t+ @ e 


(ILV¥-20b) 


yoren® = =pPReRe ee T u2eep 22m) + iajiie 
onen = —pu a as 
3P ag (xX, y,t) oie [z ( a, ( a Je 


(IV-20c) 


cs ee L and M are all complex quantities. The subscript 
a a 
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indicates the perturbation source. Garrick and Rubinow 
defined the lift force as positive down. Therefore, equating 
Eqs. (20b,c) to the negative of the corresponding components 


of Eqs. (10), one obtains, 


(| yh h 
pupae a5 : 
ket a (i) + (ed ire (L, + ar! . ole + ae 


(IV-21a) 


2 4h, h, ' 
ale oe + (Hd iS ee + oe + a, (H, + in) 


(IV-21b) 


using the definition of reduced frequency from Eq. (11-39). 





Hence, 
be - =. oe! (IV-22a) 
E = = eae) (IV-22b) 
1B = The (3 (EV-22c) 
Li = 5 att (IV-22d) 
ee == foe (IV-22e) 
M, -— Tm {t} (IV-22f) 
i — By | (I V-22g9) 
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2 
io = — 10th IV-22h 
yo ee ) 


The computer program, which is listed after the appendices 
and is explained in Section V, computes the non-dimensional 
lift force and moment defined by Eqs. (20b,c) for each blade 
of the cascade. These values are then modified by Eqs. (22) 


th 
so that the flutter analysis can be made for the n blade. 
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V. COMPUTER PROGRAM 
A. INTRODUCTION 


The finite difference equations derived in Section III-B 
and the flutter eguations of Section IV were programmed in 
FORTRAN IV and run on an IBM 360/67 computer. The computer 
program listing for the two degree of freedom case of pitch 
and plunge is given following Appendix B. 


The computational procedure is the same for pitch or 
plunge, the oniy difference between the two being the forn 
of the flow tangency equation. The two degree of freedon 
solution is simply a superposition of the individual cases 


of pitch and plunge. 


It should be noted that the subscripts of the four grid 
points of the computational molecule shown in figure 3 are 
distinct from the notation used for subscripting the 
perturbation quantities in the computer program. There, the 
Subscripts indicate whether or not a perturbation quantity 
@te oa grid point is above or below a surface, and if the 


perturbation is due to oscillation in pitch or plunge. 


In the computer program, all perturbation quantities 
have an alpha-numeric subscript (two numbers and a_e letter) 
and an index. The index identifies the grid point. The 
letter in the subscript tells whether the component is real 


(R) or imaginary (1), while the numbers indicate: 


22 (pitch) or 44 (plunge) ..... perturbation 
quantity is at a grid point which is not 
on the lower surface of a blade or wake 
slip plane. 


65 





Bout pltenleor So (plunge) -.... perturbation 
quantity is at a grid point located on 
the lower surface of a blade or siip 


plane. 
TiS , vee is the complex component of the perturbation in 


normal velocity, due to plunge, on the lower surface of a 


blade or slip piane. The subscript of the perturbation 
quantity solved for in this section merely indicates the 
real or imaginary component of that quantity at grid point 
22. 


Be. PROGRAM CONSTRAINTS 








WOO Ce, 
RESON” 





Figure 10. The compatible stagger angle. 


The computer program is quite general in its ability to 
solve a given cascade flow problem. The one geometric 
constraint that must be satisfied by the user is Eq. (Ir-d)e 
Subsonic leading edge locus of the cascade blading. The 


computer program requires that the leading edge of each 
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blade be exactly on a+ grid point. Therefore, the input 
stagger angle is modified so that STGR, shown in figure 10, 
is a non-zero integer multiple of DSTSTR aft of the initial 
Mach wave. Then the ‘compatible' stagger angle is determined 
and given in the output. The amount by which this procedure 
alters the input cascade geometry is a function of the 
gridfineness, Mach number, solidity, and distance between 
blades. See also V-G. 


The only other restrictions are on vector size. Internal 
checks are made at the beginning of a run to insure that the 
program dimensions are large enough for the input geometry 
and Mach number. If the storage area is exceeded, program 
termination occurs. Of course, the program is only valid for 


M > 1, and in the range where linearized theory is valid. 


The vector dimension reguirements for a given cascade 
and Mach number can be determineg prior to running the 
program. The vectors used for storage of all perturbation 
quantities and grid point locations along a right-running 
Mach line must be dimensioned an integer value greater than, 


1 1 
5° tan(a) {(N, q ttan(A) - Ga) Vita) + 1 (V- 1) 


1 


where, 


e *® © e grid fineness ratio 


N e « 6 e number of blades 
bid 


d e ee e distance between blades 
and @ and B have been defined previously. 
The program marches down successive right-running Mach 


lines from one grid point to the next. The location of a 


Mead point, along the current Mach line, at which the 
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perturbation quantities have just been computed is IHAVEP. 
The next grid point is labeled PHAV EDS 1. In the 
subroutines, these same two grid points are labeled J and 
I respectively. Then IHAVEP+1 = I is the index of the 
storage location within a vector array of a perturbation 


quantity or grid point location at =, or figure 3. The 
left-running Mach line which passes through grid points a 
and F is a line of constant IHAVEP = J, and the next 


21 
left-running Mach line is a line of constant IHAVEP+1 = I. 


The last grid point in the flow field is the first wake 
grid point aft of the last blade in the cascade. Along the 
left-running Mach line which passes through this grid point, 
IHAVEP = LIMIT. When the computation reaches this Mach line, 
the program Switch variables are reset, and the 
computational procedure jumps to the next right-running Mach 
line. Thus, perturbation quantity indices are from 1 to 
LIMIT+1, while IHAVEP goes from QO to LIMIT along a Mach 


line. 


The perturbation quantity vectors are single dimension 
arrays in order to reduce reguired core space. However, as 
the computational molecule requires perturbation quantities 
from two different Mach lines, it is necessary to 
temporarily switch and store the perturbation quantity 
values as scalars upon entrance to each subroutine before 


solving for ‘new' values at oe Figure 11 shows the 


Storage locations of a representative perturbation quantity 


vector as it is in the program on entrance to and exit from 
a Subroutine. 
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ENTRANCE EXIT 


HAVE | taaver-t_| 








Perturbation quantities 
from the current right- 


running Mach line. 
IHAVEP- 


| 
IHAVEP=J 


IHAVEP+¢I=I 
LIMIT +1 


Figure 11. Typical perturbation guantity 
vector storage on entrance to and exit 
from a subroutine. 















Perturbation quantities 
from the previous right- 
running Mach line. 





LIMIT +1 


The computational procedure in subroutines GENFPT, 
NEWBLD, TOP, BOTTOM and WAKE is the same. This procedure 


is; 

Step 1. The incomming perturbation guantities at 
grid point IHAVEP=J on the current Mach line Cs of figure 
3) are stored as the scalar quantities having suffixes NEW 


(UZRNEW, V4INEW, etc). Stored in the vector on the left at. 
location IHAVEP are the perturbation quantities from the 


previous Mach line eae 


Step 2. The constants K12R, K12I, etc., defined in 
Section III-B, are computed using the perturbation 
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a 


quantities at grid points Ee and De Now, perturbation 


quantities from the previous Mach line, stored at location 


IHAVEP (vector on the left) , are no longer required. 


Step 3. Perturbation quantities from the previous 
grid point of the current Mach line (U2RNEW, V4INEW, etc.) 
are stored in the vector at position IHAVEP=J (vector on 
maght) , which removes fron memory the perturbation 


quantities from the previous Mach line at grid point IHAVEP. 


Step 4. The perturbation guantities at LTHAVEPt? = I 
are computed and stored as the scalar quantities U2RNEHW, 
VULNEW, etc., and the subroutine is exited. 


This procedure is slightly different in subroutines BOTTOM, 
WAKE and NEWBLD because perturbation guantities for lower 
surface points are computed. When a lower surface point has 
been computed by one of these subroutines, the integer 
variable ICO is set equal to one. This flagged value then 
causes GENFPT, which will always be the next subroutine 
called, to use lower surface perturbation quantities in the 


computational molecule at grid point Be for the calculation 


of perturbation quantities at grid point cote 


The last grid point on all right-running Mach lines will 
lie in either the general flow field or the wake of one of 
the blades. When this grid point is reached, after returning 
from the appropriate subroutine to MAIN, the perturbation 
quantities just computed are stored. Then the program jumps 
to the next right-running Mach line. The program terminates 
flow field calculations on the first grid point in the wake 
aft of the last blade. 


Subroutine COEF stores the pressure difference across a 
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blade resulting from pitch or plunge. These values are then 
used to compute the lift and moment (about the elastic axis, 
x) for each blade. In order to be able to identify these 


pressure differences with the appropriate blade, the 


information is stored in a two-dimensional matrix. The first 
dimension identifies the blade, while the second dimension 
specifies the maximum number of grid points along the blade. 
In order to minimize the required core size, it was 
desirable not to have the first dimension equal the total 
number of blades in the cascade. Instead, a smaller matrix 
was used cyclically for all of the blades. However, this 
required that a restriction be placed on the number of 


blades that any right-running Mach line can cross. 


After the first wake point aft of a blade is reached, 
COEF2 (alternate entry point for subroutine COEF) is called, 
and the pressure difference at the end point of the blade 


(x = 1) linearly extrapolated. Then the dimensionless 


lift force and moment, as given by Eqs. (III-76), are 
computed using a trapezoidal integration technique. Once 
the lift and moment have been computed that portion of the 
matrix storing information relative to that blade may be 
reused for another blade. The maximum number of blades that 
a cight-running Mach line may cross is ICROSS, which is 
Specified in subroutine INPUT by a data statement. In order 
not to exceed the storage space allocated for these vectors 
storing pressure and moment distributions, the right-running 


Mach line which passes through the leading edge of the ({n + 


th 
meROSS + 1) blade must pass aft of the first wake grid 
' th 
point of the n blade. Also, there must be less than 


MAXPTS grid points along a blade, including the extrapolated 
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point at x A = 1. ICROSS and MAXPTS are the dimensions of 
re 


the vectors storing the pressure and moment distributions. 


Letting [| | denote the chopoed integer value of a real 


number, the restriction on blades crossed is satisfied if 


1 
ICROSS{{d[tan(@)-cot(@) VY Ax! +e} > in 2 
(V- 2) 


The right hand side of Ey. (2) gives the total number of 
grid points along the blade (which must be less than 
MAXPTS). x is the step size along a streamline (DSTISTR), 
and is egual to 


Ax = 2d[cot(a@) //e (v= 3) 


All terms were previously defined in Eg. (1). 


One final dimension that must not be exceeded is the 
total number of blades, which can not’ be greater than 
MAXBLD. MAXPTS and NAXBLD are specified by a data statement 
in INPUT in the same manner as ICROSS. Integrated lift 
forces and mcments are stored by blade number. From figure 
11, it can be seen that apparent blades are 'seen' by the 
computer. Depending on the input parameters (low Mach 
humber and high solidity), the computer may ‘see* several 
apparent blades. Therefore, MAXBLD should be at least two or 
three greater than the number of blades in the cascade. The 
required core size is insensitive to the magnitude of 
MAXBLD. 


Thus , the constraints and restrictions on program use 
eemegiven by Eq. (II-1) and Eqs. (V¥-1) - (V¥-3) and the 
comment on the maximum number of blades above. All of the 
constraints are checked internally by the program in INPUT, 


and an appropriate error termination wessage given if one or 
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more of the constraints is not net. 


C. PROGRAM FLOW AND CONTROL 


With the exception of the main program and subroutine 
COEF, all subroutines operate in a straightforward manner to 
solve the finite difference equations developed in Section 


III-B. Hence, they will be only briefly discussed. 


MAIN controls the progran flow as consecutive 
right-running Mach line are traversed. Integer switch 
variables, listed in Appendix B, are used in the the five 


primary logic tests. These are, in the order encountered, 


(1) I£ ( LCOUNT = 0 ) grid point is on the initial 
right-running Mach line 

(2) If ( THAVEP = LCOUNT ) grid point is on a blade 
or in a wake 

(3) If (IHAVEP = 0 ) grid point is on the initial 
left-running Mach line ; 

(4) If ( IHAVEP = NUMPI*NSTPTS ) and (JLINE = MAXI ) 
grid point is on the leading edge of a 
new blade 


If (2) was satisfied, then, 


(5) If ( IHAVEP > LIMW ) grid point is in the wake 


To clarify program operation and the use of these switch 


variables, the sample case shown in figure 12 is discussed. 


For M = 1.3 and the geometry given, the quantities 
listed in figure 12 are computed prior to beginning any flow 
tield calculations. ‘ ais the compatible tagger angie 


computed by the program. Subroutine INITAL calculates the 
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perturbation quantities at the leading edge of the first 


blade and initializes IHAVEP = 1, JLINE = 0 and LCOUNT = 0. ° 


MAIN then computes NEWDST and LIMIT as Shown in the figure, 


and initializes 


LIMW = NEWDST = 10 
MAXI = NGRDFN + NSTPTS = 10 
NUM = 0 
NUMOLD = 0 
NUMPI = _1 
OLDL = 0 
Ico = 0 


Now the loop which calculates the entire flow field is 


entered. 


blade 


Sit ttrrr rec rescrecri-rr-3 





INPUT 
M = 13 
€ 6 
B = 628 
d 33824 
Noid = 2 
COMPUTED 
Ax = .09365 
O = 1.356 
[2 82.7 
NSTPTS = 4 
NEWDST = 10 
LIMIT = {5 
IHAVEP = LIMIT ag 
Figure 12. Sample Cascade Problem 
IHAVEP = 0O means that the grid point is on the initial 
Jeft-running Mach line, while IHAVEP # 0 , LCOUNT = O 
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indicates that the gid point is on the initial 
right-running Mach line. In either case, subroutine MACHLN 
is called to compute the the perturbation quantities 
associated with each grid point. For computations everywhere 
except along the initialright-running Mach line, the 
temporary shuffling of perturbation guantities must te done 


prior to computing the values at the new grid point. 


After computing down the initial Mach line to IHAVEP = 
LIMIT, the program arrives at the first grid point on the 
second right-running Mach line with ITHAVEP = 0 and LCOUNT = 
JLINE = 1. All other integer switch variables are the same 
as before. At the second grid point on this Mach line (point 
1 ain figure 11), IHAVEP = LCOUNT = 1, which results ina 
call to TOP/BOTTOM (Since IHAVEP < LIMW). Prior to calling 
TOP, ICO is set equal to 1 so that the next time GENFPT is 
called (at 2) the perturbation quantities from the lower 
surface of the blade will be used in the computational 


molecule for grid point Se After return from BOTTOM, COEFi 


feeernate entry to subroutine COEF) is called to store the 
pressure difference just computed across the first blade at 


this grid point. Then IHAVEP is incremented by 1. 


For the remaining grid points along this Mach line, 
Since all of the logical tests will be false, GENFPT will be 
called by default. Again, ICO = 1 indicates that a lower 
surface point is to be used for perturbation quantities at 


Gi2d point By in the computational molecule. Prior to 
exiting GENFPT, IcO is reset equal to zero. For the 


remainder of the grid points on the Mach line. GENFPT will 
be called until LHAVEP = LIMIT. The perturbation quantities 
computed at each grid point are stored in vector locations 
THAVEP+1. When LHAVEP = LIMIT, the integer switches are 
reset to 

IHAVEP = 0 
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LCOUNT = LCOUNT+1 = 2 
JLINE = JLINE+1 = 2 
Ico = 0 
IcO is reset to zero here to allow for the case when the 
Mach line terminates on a grid point in the wake (as at 3), 


Since a call to WAKE causes ICO = 1. 


The computational procedure remains the same through the 
tenth Mach line (along A). MACHLN is called whenever IHAVEP 
= 0, and TOP/BOTTOM are called when IHAVEP = LCOUNT. GENFPT 
is called otherwise. At the completion of computation along 


(A), the integer switch variables are reset to 


IHAVEP = Q 

LCOUNT = 10 

JLINE = 10 
Ico = 0 

and the other integer variables are as before, 

LIMNW = 10 

NSTPTS = 4 
NUMPI = 1 
MAXI = 10 


The computational procedure down Mach line (B) is the same 
until the grid point at (4), where, 

IHAVEP = (NUHNPI) (NSTPTS) = 4 
and 

JLINE = MAXI = 10 

Satisfaction of this logical test envokes a call to NEWBLD 
to compute the perturbation quantities at the leading edge 
of the second blade. NEWBLD also also increments NUM so that 
NUM = 1 on return . Next, COEF1 is called to store the 
pressure difference at the first grid point on the second 
blade. Next,the integer switches are reset to, 


ICO = 
EeOuNT = ITHAVEP + NGRDFN = 10 


LI = 
IHAVEP = IHAVEP + 1= 5 
OLDL = IHAVEP = 5 


NUMPT NUM + 1 = 2 
MAXI = NGRDFN + (NUMPTI) (NSTPTS) = 14 
NUM = NUM ~ 1 = 0 


Computation down (B) continues until the grid point at (5), 
with GENFPT keing called for all intervening grid points. 


For the - first grid point after (4), IcO = i, and was reset 
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equal to zero in GENFPT. At (5), IHAVE = LCOUNT = 10, and 
TOP/BOTTOM are called. COEF1 is called on return fron 
BOTTOM, and GENFPPT is called for the remaining grid points 
until IHAVEP = LIMIT. Then the integer switch variables are 
reset to 


IHAVEP 


= 0 
NUM = aad = 


1 
ICO = 
JLINE = JLINE + 1 = 5 
Since more than one blade has now been encountered, NUMOLD > 


0, and it iS necessary to reset the following integer 
switches so that the grid point on blade two can be 


identified. 
LCOUNT = OLDL = 5 
LIMW = (NUMOLD) (NSTPTS) + NEWDST = 14 
OLDL = OLDL + 1 = 6 
The computation now starts over at the initial 


left-running Mach line (Mach line C), where MACHLN is called 
for IHAVEP 
where IHAVEP = LCOUNT = 5, GENFPT is called. At (6), ICO is 
. set egual to one and TOP/BOTTOM are called. Next, COEF1 is 


called to store the pressure difference just computed. Upon 


Ge Thereafter, untii the grid point at (6) 


return from COEF1, since NUM > 0O, the integer switch 
variables must be reset so the grid point on the upper 
surface of the next blade (or wake, in this case) can be 
identified. Thus, 
NUM = NUM - 1 0 
LCOUNT + NGRDFN 11 
LIMW = LIMW ~- NSTPTS = 10 


GENFPT is then called to compute the perturbation quantities 


at grid points until the program reaches (7), where IHAVEP = 
LCOUNT = 11. This time, IHAVEP > LIMW , which means that 
WAKE should Le called instead of TOP/BOTTOM. Again, ICO is 
set egual to one before calling WAKE. After return from 
WAKE, COEF2 is called. As this is the first time that it 
has been called for this blade (i.e., this is the first grid 
point in the wake aft of blade one), the intecration of the 


pressure distribution for the lift and moment is done. 


Wi 





Since this is the first call to COEF2 aft of blade one, 
IWRITE =. 1 on return to the main program and the pressure 
distribution for this blade can be output, if desired. Next, 
a check is made to determine whether or not this was the 
last grid point in the flow field, or if IHAVEP = LIMIT. If 
not, the calculation loop is reentered, with ICO = 1. GENFPT 
is called to compute the remaining grid points until IHAVEP 
= LIMIT, where the integer switches are reset to 


IHAVEP = 0 
NUM = NUMOLD = 1 


Ico = 0 
JLINE = JLINE + 1 = 6 
and as NUMOLD > 0, 
LCOUNT = OLDL = 6 


LIMW = (NUMOLD) (NSTPTS) + NEWDST = 14 
OLDL = OLDL + 1 = 7 


The computational procedure down Mach line (D) is the 
same as for (C). The call to COEF2 at (8) will result in an 
immediate return to the main program with IWRITE = 0 as this 
will be the second time (more accurately, not the first 
time) that COEF2 has been called for blade one. The 
computation proceeds in the same manner, down successive 
right-running Mach line, until reaching point (9) on Mach 
line (E), where, 

ITHAVEP = (NUMPI) (NSTPTS) = 8 

and 

JLINE = MAXI = 14 
which initiates a call to NEWBLD for the assumed blade. 
Termination of the flow field calculations is indicated when 

NUM > NUMBLD - 1 = 1 
and 

IHAVEP = LIMIT 

at point (10). As (10) is also the first wake grid point aft 
of blade 2, COEF2 is called to compute the lift and moment 
for the second blade before terminating the flow field 


calculations. 


78 





D. SUBROUTINE COEF 


Subroutine COEF is a multiple entry subroutine used to 
store the pressure difference distributions for each blade 
as it is encountered, to extrapolate for the pressure 
difference at the trailing edge of each blade, and to use 
these distributions to compute the lift and moment due to 
pitch and plunge. If the pressure distribution fora 
particular blade is required, output is obtained by issuing 
a write statement after return from COEF2 in the main 


progran. 


COEF is the first section of this subroutine. Called 
after return from INITAL, it initializes all arrays used in 
the subroutine and stores the pressure difference for the 
first grid point on the first blade. The pressure 
difference is then used to compute the incremental moment 


about the elastic axis. 


COEF1 is called after every call to TOP/BOTTOM or 
NEWBLD, and stores the pressure difference with the real 
component referenced to the blade initial position, and the 
imaginary component referenced to the mean position, blade 
pitching leading edge up. Refer to _ Eqs. (E1I=— 70) = 
(III-75). The pressure differences are then used to compute 
the incremental moment resulting from each discrete pressure 


difference. 


COEF2 is called after every call to WAKE. The first 
time this section is called, it initiates the calculation of 


th 
the lift force and moment for the an blade, after 


extrapolating for the pressure difference at the trailing 


edge of the blade. At the same time, the integer switch 
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ID(n) is set egual to one. Subsequent calls to COEF2 for the 
th 


n blade wake result in an immediate return to the main 
program, after setting IWRITE = 0. IWRITE is an integer 
switch variable used to indicate when the pressure 


‘distribution for the entire blade is complete. For a given 
blade, IWRITE = 1 only after the first call to COEF2, where 
the pressure difference at the trailing edge of the blade is 
determined. At all other times, IWRITE = 0. After the lift 


force and moment have been computed, the section of the 


matrix used EOL storage of the n blade pressure 

distribution can be reused for another blade. Upon return 
th 

to MAIN, the pressure distribution for the n blade can be 


output, if desired, if IWRITE = 1. 


E. OTHER SUBROUTINES 


The remaining subroutines used to calculate perturbation 
quantities in the flow field are TOP, BOTTOM, GENFPT, 
NEWBLD, WAKE, and MACHLN, all called by the main progran. 
Each subroutine solves the applicable finite difference 
equations derived in section III-B. Subroutine WAKE uses IBM 
subroutine SIMQ, a simultaneous equation solving routine, to 


compute the perturbation quantities at F The solution 


22 
procedure in each subroutine is the four step process 
described earlier in this section, and it is carried out 
first for pitch, and then for plunge in each subroutine. 
Each subroutine defines several constants, a few of whict 


are common to both solutions. 
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Subroutine INPUT reads the input data, computes the 
compatible stagger angle, and the geometric cascade 
parameters used by other subroutines. Here the input data is 


checked to insure that available core size is not exceeded. 


Subroutine INITAL uses the equations derived for the 
initial Mach lines to compute the perturbation quantities at 
the leading edge of the first 
blade. In addition, JLINE, LCOUNT 
and IHAVEP are initialized. 

DELTAS 

Finally, subroutine COMPXY 
computes the grid point DSTSTR 
coordinates in the characteristic 


Figure {!3. Characteristic 


mesh in terms of the the geometry mesh geometry 


computed in INPUT and shown in 

figure 13. The calculation procedure for points on _ the 

right-running Mach line is 
X (I) 


X(I) + HDSTRL 


¥ (I) Y(I) -— TRNGLH 


The first grid point on a right-running Mach line is 
Ketek (1) + ADSTRL 


Y¥(1) = ¥(1) + TRNGLH 


Actually, the only coordinate used by the program is X(I). 


Subroutine FLUTER is called for two degree of freedon 
flutter solutions when the input variable IFLUT = 1, and 
uses the equations derived in Section IV, where the _ single 
degree of freedom equations were given also. Flutter 
investigation is quite expensive in terms of computer time, 
especially if the program is designed to iterate to an 
actual flutter speed, subject to some error tolerance. To 
Pepeamize the 'first blade effect", the flutter calculations 


must be conducted for a relatively high numbered biade where 


81 





the magnitude of the lift forces and moments have nore 
nearly approached an assymptotic value. This must be done 
because an actual turbomachine does not have a 'first' blade 
and the aerodynamic forces on the blades are assumed _ to be 
periodic. To determine a critical flutter speed for a given 
cascade geometry and Mach number, the minimum flutter speed 
must be determined, which requires an iteration on the 
interblade phase angle as well. Since each iteration 
through the flutter routine requires a recomputation of the 
entire flow field, it is apparent that computational time 
mounts rapidly. 


In an attempt to reduce the amount of computer tine 
required for the flutter investigation, the conputer program 
is set up to compute and output the roots of the flutter 
determinant versus the reciprocal of the reduced frequency 
(1/k). These results are then plotted by hand, using the 


Theodorsen method to find a flutter solution. 


€1/k) is incremented by an amount STEP, an input 
Variable. If, for a given value of (1/k), the real or 
imaginary component of the solution is negative the program 
will terminate, as the dimensionless flutter frequency would 
be imaginary. However, if the discriminant of the quadratic 
for Re {X} is negative, then (1/k) is first decreased by 
STEP, and then increased by one-fifth that amount. The flow 
field is then recomputed, and FLUTER called again. Note that 
this procedure assumes that the reduced frequency used for 
the previous solution produced a positive real root. As the 
real component of the flutter determinant is quadratic, a 
negative discriminant indicates that the value of (1/k) is 
past the vertex of the parabolic curve of /Re {Xx} in the 
solution plane (see figure 9, Section IV). Each time 
through the flutter subroutine, the square root of the three 
Peot= Of the flutter determinant and the cortysponding valu 


of (1/K) are output. Because of the manner in which the 
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case of a negative real root is handled, and also because 
the minimum flutter speed occurs at the smallest value of 
(1/k), the program should initially be given a reduced 
frequency on the order of three or greater. In all cases 
investigated in this report, all three of the flutter 
determinant roots were positive for reduced frequencies of 


that magnitude. 


F. INPUT AND OUTPUT 


Subroutine INPUT reads all input data from two Namelist 
input records, NAM1 and FLTR. The input variables for NAM1 
are, 

NGRDFN - the grid fineness ratio (e). The number 
of increments along a Mach line between 
two blades. 

FSTRMN - the freestream Mach number. 


i 


REDFROQ - the reduced frequency, defined in Eq. 


(OIE BS) 

XSUBO - the elastic axis position (x) 

TNWDST - vertical distance between adjacent 
bliades. 


STGANG - the input stagger angle (8). The 
compatible stagger is computed in INPUT. 


FAZE - the interblade phase angle (in degrees). 
NUMBLD - the nunber of blades ee. rg the 
cascade. 


IPRES1 - number of the blade user desires specific 
information about, 1.@., pressure 
Gastrabution, Litt force and moment 
values, etc. 


Same as TPRES 1. 


L 


HPRES2 


Variables in Namelist FLTR are all used in the flutter 


routine. They are, 


RMUU - the blade density parameter (p). 
XSUBA - location of the center of gravity 
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measured from the elastic axis, positive 


aft (Xx ,)- 
RSUBA - the blade radius of gyration (r,)- 
WHWA - ratio of natural frequency in plunge to 

@w 
the natural frequency in pitch () 
a 

STEP - increment for (1/k). 
IFLUT - set egual to one (iFLUT=1) if a flutter 


solution is desired. 


The input data for the sample program discussed in part 


C of this section would be: 


MHD OOH 
leh tlcketuolssiaclicl h- 
Db 
Te 


20 Bsa cach aymaxs 
Var BOCHUM hro 
iH 
Sy) 


Each card begins in column 2, and the comma follows 
immediately after the last digit. Spaces between the last 
digit and the comma are interpreted by the computer as 


zeros. 


Output can be modified at the users discretion. With the 
program as listed after Appendix B, pressure distributions, 
dift forces and moment and the magnitude and phase of the 
complex lift and moment for each blade can be output. In 
addition, if appropriate ‘write' statements are put in each 
of the subroutines, the entire flow field may be printed. 
This option produces sizeable quantities of output, and is 


useful primarily for program debugging. fic output Eros 
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flutter solution runs was described in Section V-D. 


G6. GENERAL 


Because the program computes a compatible stagger angle, 
the cascade geometry actually used in the computer progran 
may be slightly different than that input by the user. For 
large values of the grid fineness ratio,e , the differences 
are insignificant. However, it is normally preferable to use 
a coarser grid size to minimize computer time. As a result, 
if a specific cascade is to be input, then an iteration on ¢ 
Should be done to find the smallest value of ¢ which 
produces the compatible stagger angle that is in closest 
agreement with the input stagger angle. If this procedure is 
not followed, the cascade computed for will be slightly 
different than desired, and therefore, the reflection points 
on the blades will be different than expected. 


The computer program requires approximately 175K bytes 
of storage (G compiler) or 135K bytes of storage (H 
compiler) when run on an IBM 360/67 computer. Compile time 
is approximately twenty-five seconds. Actual CPU time per 
Cun is dependent on the input geometry, grid fineness and 
number of blades. Sample times are given for each of the 


results discussed in Section VI. 


Initial progran checkout and verification can be 
accomplished in several ways. By setting the reduced 
frequency equal to zero, the well known result for the lift 
on a flat plate at unit angle of attack is obtained (for the 
first blade). Also, tabulated results appearing in the 
report by Garrick and Rubinow [Ref. 9] for single blades 
require less than five seconds of computer time to generate 


using a grid fineness of 100 or less. 
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VI. RESULTS 


Initial work was done for single degree of freedon 
torsion. All cases discussed in this section were run with 


zero Structural damping (g = 0). Pressure distributions 
and flutter solutions were computed and compared with Snyder 


(Ref. 19}, who used a finite difference procedure developed 
by Verdon [Ref. 23]. Brix [Ref. 2] showed that considerable 
discrepancy existed between the pressure distributions given 
by Snyder and those computed by the method of 
characteristics, especially in the region of reflection 
points. As Stated by Verdon [Ref. 25, p24], 

The finite difference marching procedure used .. . does 

mot Tecognize pressure discontinuities as such, but only 

approximates them by sharp changes in the _pressure 

distribution. This numerical approximation also gives 

Lise to the ripples which opneet ay ele ae cs bata 

o 


ressure |. distributions ownstreanm the 
iscontinuities. 


The flutter comparisons shown in figures 13 and 14 are 


for cascades C, D and E, described below. 


Cascade C Cascade D Cascade E 

B = 67 B = 60 j= 70 
M € Tine € ‘Time € Time 
1.3 24 70 55 75 26 85 
1.4 27 50 39 60 20 45 
LAS) 27 45 LQ 50 22 45 
eG 27 35 31 DS 20 25 
1. 7 28 30 Bal 4Q 22 30 
1.8 28 25 38 20 24 30 


The entry in the Time column is the actual CPU time, in 
seconds, required for one iteration through the flow field 
for a ten blade cascade. Figures begin on page 91. 


Meieesthree —cGascades haye a solidity of 1.33 and a 
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mid-chord pitch axis (x 4)» The grid fineness, e€«, was varied 


in order to be able to match the geometry of each cascade as 
closely as possible (refer to V-G), while maintaining at 
least twenty-five discrete grid points along each blade. 


Both results assumed zero structural damping (gi). In 
addition, the method of characteristics solution used an 


assumed value of 500 for the blade density parameter p, and 
a radius of gyration, Tye of 0.5. The value for p» was felt 


to be reasonable for modern turbomachinery blading [Ref. 5]. 


These figures, 14 and 15 show the reduced frequency which 


occurs at the minimun flutter speed, k earand the 
Cra 


interblade phase angle at that point, 8 —, versus Mach 
(yealie 


number. Method of characteristics results are based on the 


tenth blade. Agreement is not exact. However, in light of 
the previous comment on the computational scheme used by 


Snyder, the agreement is felt to be reasonable. 


In figures 16 to 20, Similar results are shown for 
Cascade C with elastic axis aS a parameter. Again, p= 500, 
and r = 0.5. 

a 


Figures 21 to 24 present flutter stability boundaries 
for Cascade C. These are plots of elastic axis position 


versus Mach number for which oe = 0. The shaded portion of 
each graph is unstable, i.e., ar < 0. The analysis was made 
for 30° increments in the interblade phase angle, §. Figures 


21 and 22 showed instabilities at only the values of 
Shown. In figure 23, the range of interblade phase angles 
over which instabilities existed increased, aud in figure 


24, instabilities occurred for all values of examined. 
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Thus, these four figures show that as the reduced frequency 
is increased (1/k decreased), cascade operation becomes nore 
stable. 


When conducting flutter investigations, it is impcrtant 
to find the highest value of reduced frequency for which the 
cascade becomes unstable, as this will produce the lowest 
flutter speed for a given geometry and interblade phase 
angle. Figure 25 shows the variation of a With (1/k). As 


can be seen, if the step size in (1/k) 1s too large, the 
first point at which ar < 0 could be missed (note inset in 
the figure). As a result, the reduced frequency at which M 


= 0 at B would result in an incorrect (too high) flutter 


speed. Riso, figure 25 shows the condition where an just 
crosses the axis. When the phase angle was’ slightly 
decreased, the curve of Ze Was not quite tangent to the axis 
at A, but crossed the axis near B. Figure 26 shows the 


resulting jump in the flutter speed and reduced frequency at 
putter. 


A limited amount of investigation was conducted to 
determine the dependence of flutter speed on blade number. 


Although S is quite oscillatory for the first few blades 
and does not approach an asymptotic value until after the 


tenth or fifteenth blade, preliminary results, based on a 
limited number of cases, indicated that the flutter speed 
and frequency were fairly constant. An example is shown in 
figure 27 for Cascade D with M = 1.7, 8 = ~45 , and ¢€ = 25. 


Comparison with results by Verdon [Ref. 25] for an 


infinite cascade wasS made using the two de.vee of freedom 
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computer program listed at the end of this report. The 
cascades A and B referred to in the figures are shown and 
described in figure 28. Pitch results are always for the 
fifteenth blade, and plunge results, for the fourteenth. In 
figures 29 through 49, the time dependence is in terms of 
the compressible reduced frequency as defined by Verdon, 


where k = xe Pe Corresponding to Verdon's notation, 
(3 -_ 


Im{M} > O implies unstable operation for the pitch case. 
Plunge oscillation was always found to be stable. [In 


addition, figure 29 shows convergence of the norgm of M 
erm = {js | = M eM ) versus blade number for four 
a 


different values of 86. These phase angles are all in the 


region for which Verdon was unable to obtain a_ solution. 
The curves are shown aS being continuous for clarity of 
presentation, but results were only obtained at discrete 
integer values of blade number, of course. Each curve 


reguired approximately nineteen minutes of computer time. 


Figures 30 through 33 are plots of the real and 
imaginary component of the pressure distribution acting on 
the upper and lower surfaces of the blades for the two 
cascades. Figures 34 through 37 show the pressure difference 
across the blade (lower - upper) for two different 


interblade phase angles. Agreement was excellent. 


Figure 38 defines the symbols used in figures 39 through 
49, where lift and moment loops are Obtained by plotting the 
real and imaginary components of the 1ift force and the 
moment with the interblade phase angle aS a parameter. The 
method of characteristics was able to calculate solutions 


throughout the range of interblade phase angles, 0 < & < 
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360. Again, excellent agreement was obtained. 


In figures 30 through 49, results for cascade A required 
48-50 seconds of CPU time per fifteenth blade solution, and 
for cascade B, 100-110 seconcs per fourteenth blade 
solution. 


The last figure, 50, shows a flutter solution for 


cascade A assuming the following parameters: p = 500, tat 
=n 
0.5, x = 0, and — = .707. 
a on 
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& Snyder and Commerford 
(Ref. 19) 
C Present study (Nig = 10) 
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gure 15. Interblade phase angle at minimum flutter speed 
versus Mach number (2 =500,r, =0.5 ); 
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Figure 16. Reduced frequency (kcrit) and interblade phase 
angle ®crit) at minimum flutter speed versus Mach number 
005m ao. Ji: 
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angle (crit) at minimum flutter speed versus Mach number 
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Figure 18. Reduced frequency (k,,4;) and interblade phase 
angle (Sor;¢) at minimum flutter speed versu. Mach number 


( #=500,r,=0.5 ). 
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Figure 19. Reduced frequency a) and interblade phase 


angle (8.;;;) at minimum flutter speed versus Mach number 
( H=500,r,=0.5 ).. 
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Figure 20. Reduced frequency (keri¢) and interblade phase 


angle (one) at minimum flutter speed versus ‘ach number 
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Flutter boundaries versus elastic axis position 


Figure 22. 


(shaded portion is unstable). 
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Figure 23. Flutter boundaries versus elastic axis position 
(shaded portion is unstable). 
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Figure 24. Flutter boundaries versus elastic axis position 
(shaded portion is unstable). 
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Figure 25. 


Oscillation of the imaginary component of the 
moment due to pitch (My) with reduced frequency. 
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Figure 26. Variation of flutter speed and reduced frequency 
at flutter with interblade phase angle ( 4=500,r,=0.5 ee 
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Figure 27. Non-dimensional flutter speed and reduced frequency 
at flutter versus blade index ( 2=500,r,=0.5 ye 
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Figure 29. Convergence of the norm of the moment due to pitch 
oscillation with blade index as a function of the interblade 
phase angle. 
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Figure 30. Comparison of the pressure distribution resulting from 
pitch oscillation of the fifteenth blade with the infinite cascade 
results of Verdon (Cascade A). 
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Figure 31. Comparison of the pressure distribution resulting from 
plunge oscillation of the fourteenth blade with the infinite cascade 


results of Verdon (Cascade A). 
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Figure 32. Comparison of the pressure distribution resulting from 


pitch oscillation of the fifteenth blade with the infinite cascade 
results of Verdon (Cascade B). 
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Figure 33. Comparison of the pressure distribution resulting from 
plunge oscillation of the fourteenth blade with the infinite cascade 
results of Verdon (Cascade B). 


110 





Cascade A 
Pitch Xer= 0/0 
ke = LO 
—— Verdon (Ref.25) 


O ® Present study 





a: 
5:0 
=> , fia ae ian ie : i 
face ge 
£ _ 5 
-5 Sa 
6) 2 4 6 8 1.0 
X 


Figure 34. Comparison of the pressure difference distribution (lower - 
upper) resulting from pitch oscillation for the fifteenth blade with 
Pemintinite cascade results of Verdon (Cascade A). 
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Figure 35. Comparison of the pressure difference distribution (lower - 
upper) resulting from plunge oscillation of the fourteenth blade 
with the infinite cascade results of Verdon (Cascade i) < 
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Figure 36. Comparison of the pressure difference distribution (lower - 
upper) resulting from pitch oscillation for the fifteenth blade with 
the infinite cascade results of Verdon (Cascade B). 
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Figure 37. Comparison of the pressure difference distribution (lower - 
upper) resulting from plunge oscillation of the fourteenth blade 
with the infinite cascade results of Verdon (Cascade B). 
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EXPLANATION of SYMBOLS 


Method of 
Characteristics 8 Verdon 
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& 90 & 
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e 8 +30 


The solid curve connecis points determined by 
Verdon (Ref. 25). The deashed curve connects 
points determined in the present study. For 
regions of S&S with common resulfs, the solid 
curve only is shown. 


Figure 38. Definition of symbols used in figures 39 through 49. 
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Figure 39. Comparison of the moment coefficients due to pitch 
oscillation for the fifteenth blade with the incinite cascade 
results of Verdon (Cascade A). 
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Figure 40. Comparison of the moment coefficients due to pitch 
oscillation for the fifteenth blade with the infinite cascade 
results of Verdon (Cascade A). 
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Figure 41. Comparison of the moment coefficients due to pitch 
oscillation for the fifteenth blade with the infinite cascade 


results of Verdon (Cascade B). 


118 





Cascade 
Pitch o = OS 
—— Verdon (Ref. 25 ) 
Present study 


unstable 


stable 





“1.0 we, .O 5) 1.0 
Re [ Ma] 
Figure 42. Comparison of the moment coefficient: due to pitch 


oscillation for the fifteenth blade with the infinite cascade 
results of Verdon (Cascade B). 
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Figure 43. Comparison of the lift coefficients due to plunge 


oscillation for the fourteenth blade with the infinite cascade 
results of Verdon (Cascade A). 
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Figure 44. Comparison of the lift coefficients due to plunge 
oscillation for the fourteenth blade with the infinite cascade 
results of Verdon (Cascade . A). 
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Figure 45. Comparison of the lift coefficients due to plunge 
oscillation for the fourteenth blade with the infinite cascade 
results of Verdon (Cascade A). 
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Figure 46. Comparison of the lift coefficients due to plunge 
oscillation for the fourteenth blade with the infinite cascade 
results of Verdon (Cascade B). 
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Figure 47. Comparison of the lift coefficients due to plunge 
oscillation for the fourteenth blade with the infinite cascade 
results of Verdon (Cascade 8B). 
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Figure 48. Comparison of the lift coefficients due to plunge 
oscillation for the fourteenth blade with the infinite cascade 
results of Verdon (Cascade B). 
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Figure 49. Comparison of the lift coefficients due to plunge 
oscillation for the fourteenth blade with the infinite cascade 
results of Verdon (Cascade B). 
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Figure 50. Non-dimensional flutter speed and reduced frequency at 
flutter versus interblade phase angle for Cascade A ( j2=500,7r,= 
On 1). 
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VII. CONCLUSIONS 


A method of characteristics procedure has been developed 
for the solution of the complete linearized unsteady 
aerodynamics of a flat plate cascade having a subsonic 
leading edge locus in a_ supersonic flow field. Pressure 
distributions were computed and the resuiting lift forces 
and moments used to investigate flutter boundaries and 
stability limits for a two-dimensional finite cascade, which 


can have arbitrary stagger angle and interblade phase angle. 


Agreement with the infinite cascade results by Verdon 
were excellent. However, further verification and comparison 
with new methods being developed by others [Refs. 12, 20,26 ] 
should be conducted. Of primary importance would be the 
development of a periodic solution applicable for the 


infinite cascade. 
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VII. CYLINDRICAL SHELL 


This developement is closely based on a previous paper 


by Platzer, Brix and Webster [Ref. 17]. 
A. PROBLEM FORMULATION 


Consider the supersonic flow of an inviscid, adiabatic 
gas through a circular cylinder whose axis is aligned with 
the freestream, and whose surface is rigid upstream of the 
origin. The cylinder geometry is shown in Piguie ss i. 








Figure 51. Cylinder Geometry 


Downstream of the origin, a-flexible panel of length L 
is assumed to undergo small amplitude vibrations. The 
governing equation for this problem is the Linearized, 


unsteady velocity potential equation in eylinudrical 
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coordinates, $(X,Rk,9,T), 


1 1 M 1 
1-He + oS + — - 2- -—— = 0 
nee P2R R°R Re 99 oxt ce’ rr 
(VIII- 1) 
The linearized form of the pressure coefficient is, 
C = é p d Vit i— 2 
p Die T u, x ( ) 


The linearized boundary condition on the cylinder inner 
surface is 
dh dh 
Sam } ih 
*ale=r, a2” “odx 





, O-cey < 1 (VIII- 3) 


The boundary condition at r = 0 iis discussed separately. 


foenq. (3), h = h(X,8,7) is the equation for the panel 
surface. One other condition that must be satisfied is the 
Sommerfeld radiation condition; i.e., waves must propagate 
away from the source of disturbance. For a derivation of the 
above three equations, refer to Bisplinghoff, Ashley, and 
Halifman, “Aeroelasticity" [Ref. 1]. 


To non-dimensionalize the above equations, lengths are 
referred to flexible cylinder length and time to flexible 


cylinder length divided by the freestream velocity. Then, 


1 1 
(1-M2) 6 + +-—p +—O —- 226 - M2 =) 
XX ese je se re 69d brats tt 
(VIII- 4) 
ee 76 + } (ates 2) 
p te x 


subject to, 
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So = 0 <.x < 1 VIII- 6 
“ c=R ot ax * a = ( ) 


The non-dimensional shell deflection is assumed to be of the 


forn, 
= ikt 
h(x,9,t) = Z(x)cos(nd@)e (VIII- 7) 
where, 
wh 
k =— (VIII- 8) 
us 
and n is the circumferential mode number. Then the 


perturbation potential must be of the forn, 


ikt 
o(x,c,8,t) = #(x,r)cos(n9)e (Viit- 29) 


and the resulting equations are, 


1 né@ 
(1-4) ¢ + g Ome ol cae) Oo = 
Xx EE oa re 
aun = (0 (VIII-10a) 
x 
Gn= = 2g + ike) . (VIL I-10b) 
p X 
C is the pressure coefficient amplitude. The boundary 
P 
conditions are, 
Seer kz ORK x <1 (VIII-11a) 


r EHS x 


The boundary condition at the axis is dependent on the 


circumferential mode number on. Thus, 


even (VIITI-11b) 


i 
© 
« 
th 
O° 
La 
S) 


odd (VIII-11c) 


i 
i) 
2 
th 
fe) 
Le 
ro] 
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Eqs. (10) and (11) are the basic equations which will be 
used in the cylindrical shell analysis. 


B. METHOD OF CHARACTERISTICS 


For M > 1, Eq. (10) is hyperbolic and has real 
characteristics which satisfy the ordinary differential 


equation 
(M2-1)d2r - d2x = 0 (VIII- 12) 


For ds the differential arc length along a _ characteristic, 


d@s = M@der, and, 





ax JM2-1 

eae (VIII-13a) 
ds M 

ar 1 

a +t (VIII-13b) 
ds M 


The derivative along a characteristic of an arbitrary 


munetion F(x,r) is, 


aF OFdx OFdr 





— 5 eed (VIII- 14) 
ds Oxds drds 
Then, 
aF /M2-1 1 
=a r 4 =F (VIII-15a) 
- ds, Mi x Mor 
dF Me-1 1 
Po=— =*-—-P -—F (VIII-15b) 
@ ads M x Mr 


where ds and ds. are the differential are lengths along 


the left and right characteristics respectively. The second 


(cross) derivatives are, 
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1 
=F = =f (l-e) F + F = =f VELI= “Us 
ep ate 2 Ee “ ( ) 


Brom Bg. (15), 


M 
P= os oe VitI-17a 
AUER ee 
ES = : F F V 17b 
i 2) (VII I-17b) 


Now, with Eqs. (16) and (17), the basic equations are 








written, 
1 : né@ 
Ce ees 5ae Beater at {ke = SEE )6 
) Ke 
~ a aeea 6. - Be ) (VIII-18a) 
a g) Z ikZ VIII-18b 
= — ~ = + 1 =_ 
a 5 oe : a ( ) 
C =- 2g + ikg] (VILI-18c) 
Pp x 


C. FINITE DIFFERENCE EQUATIONS 


The finite difference equation for the points P and § 


along the right~running Mach line is, 








_ As a né 
BS) — 6L(P) = erica Stole io ks = Ee )g 
@ _ KMAs 
S "Re | ae = B, ) (Vili- 19) 
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* 


Note that in this and subsequent equations, ds | = ds, =As. 


For points Q and S along the left-running Mach line, 








* As ; : né 
B45) ~ eg WO) sores = Bit As{ | ea es )¢ 
_ kMAs 
~ i ee | B Ae ) (VIII- 20) 


The computational molecule is shown in figure 52. 





Figure 52.,. Characteristic grid with 
computational molecule. 


The values for @, an and B, are averages between 
adjacent grid points along the respective Mach lines, so for 


the left-running Mach line, 


# 


1 
monet 2, (2) * 2, (5) ] TE 


g 


; 
Mateo) 6) 1 URES ANS) 


and ¢(S) is obtained by integrating along the Mach line. 
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a 


1 
#(S) = 2(0) + 59, (Q) + ¥, (3) Jas (VIII~21c) 


1 
ld Ble) + $(Q) ] (VITI-214) 


Similarily, along the right-running Mach line, 


1 
= = per, te) + B(S)] (ET i227) 
[i pies y B (5) ] (VIII-22b) 
1 
£(S) = A(P) + Beat) + #_(S) JAs (VII I=22c) 
1 
2B ae UZie) + g(S) J (Vitis 22d) 


Eqs. (21) and (22) are substituted into Eqs. (19) and (20). 


The result is a system of two equations for % (5) and B,(S) 


of the forn, 


i) 
O 


Ag (8S) + BB, (S) (VIII-23a) 


iT} 
ry 


Dg (S) + EB, (S) (VIII-23b) 


Eg. (23a) applies to the right-running Mach line, where 


AS . kMAs 
+ 











A = = 2 eared VIII-24a) 
4hMr (S) 2/u2-1 
3 As ( Ke ne ee _ kMAs 
= = ee —— t+ 1 
4Mr (S) M2r2 (S) Z 2/ 2-1 
(VIII-24b) 








conn, Oe ee ee ee 
2 4Mr(P) M2r2(S) 2 2/m2-7 
+g (PL 14 AS ies 
1 4Mr (P) 2/M2=-1 
n@ né@ As 
+ 2 (P) (ke 5= W252 (2) ea ke = ner2(s) ? in 
(VIII~24c) 


Eq. (23b) applies to the left-running Mach line, with, 





— As oe ne ee eee 
4Mr (S) M2r2(S) 2 2/i2e=1 
(VIII-25a) 
As _ kMAs 
= pee iTS (VIII-25b) 
P= 6 (0 - pom - (ke + 1 
1 4Mr (Q) M2r2(S) 2 2/M2=1 
as _ kMAS 
+2 Oy 1 = iene) rar 
- SQL ( cya ) + ( ke pea Tes 
M2r2 (Q) M2r2(S) 2 
(VIII-25c) 


Egs. (23) are used to compute g and its derivatives at a 


general flow field point. 


Along the initial right-running Mach line, ¢g = #, = 


0, and Eq. (23a) reduces to, 


Bg, (S) = C/A (VITI-26a) 
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and C simplifies to, 


XS _ kMAs 


© 39 (GN 0 ae 


—— VIII-26b 
4Mr (P) 2/M2=1 ( 


3. Inner Surface Of Outer Cylinder 


ne re ee a ee ea eee 


The flow tangency condition prescribes the normal 


velocity at r= Ro: Thus, 
| e 
g =-—{g (S) - @_(S)]=2 + ikz (VII1I-27a) 
r}r=R, 24 2 x 
For the initial point on the panel, B, = 0, and 
fs 
g@ (8S) =-[Z + ikZ} (VIII~-27b) 
1 M x x=0 


Thereafter, the equations for the left-running Mach lines 


are used to solve for @, an and B.: Hence, 


2 
2 (8) = (es oe + ikZ)EJ/(D + B) (VIII-28a) 
x 
2 
A) = ips (ze + ikZ)DJ/(D + B) (VIII-28b) 
Se = 26 + ikg] (VIII-28c) 
Pp x 


Then, ce and B are used to compute g andc. 
x Je 


4. Boundary Condition at the Axis 


The boundary condition at r = 0 was approximated by 
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assuming a very thin solid cylinder with radius r = R. along 
i 


the axis, extending the length of the cylinder. The form of 


the boundary condition applied at r = R. is dependent on the 
a 
circumferential mode number n. For even values of nh, 


4) = 0 nh even (VIITI-29a) 


6 (S) = C/(A + B) = (5) (VIII-29b) 


For odd values of n, the prescribed boundary condition is 

Cc = 0. For the steady case (k = 0), this condition reduces 

to ¢ = Q. Then applying the boundary condition with Eqs. 
x 


(23a) and (17), 


ic 0 = g (S) + ikg(S) n odd (VITI-30a) 
p c=R, x 
1 
O(S) = B(P) + Bele ABU + (5) JAs (VIITI-30b) 
i , 
Bete! = ica + $65) ] = - 1ikg(S) (VITI-30c) 
Ag (5) + ea) =C (ViTI—-30d) 


Defining, 


5 : ie VIII-31a) 
Zvi-— | si 2 : 
As 
aR (2) + (36, (2) ] (VILI-31b) 





Then, 


M OB 
2/M2-18 J 


(Dea ee Bei (VILI-31c) 


g, (S) [ ¢ - ag (S)y YB (VITI- 314) 


Now $(S), @ (S), and g (S) are computed using Eqs. (17) and 
x je 


(222). 


The initial right-running Mach line is assumed not 
to reflect from the assumed inner cylinder. 


D. GENERALIZED FORCES AND THE PRESSURE COEFFICIENT 


Having oktained the pressure distribution over the inner 


surface of the outer cylinder using the pressure coefficient 


amplitude, the generalized aerodynamic force Q is found 
mr 
fron, 
; i 
Ome eC (x) (x) dx (VIII-33a) 
nr 2 pm 1G 
oO 
where, 
Z(x) = aw (x) (VIII-33b) 
13 


and C (x) is the pressure distribution resulting from the 
po 


th th 
ri axial mode deflection and W (x) is the oc axial 
xr é 


deflection mode. Thus, one would be the generalized force 


due to an axial deflection of at least two modes of the 


ig 





forn, 
ZX). = AYO) + ALY) (VIII-34a) 
and, 


C  (x)¥ (x) ax (VIII-34b) 
Pi 2 


The above eguations were programmed in Fortran IV and 
run on an IBM 360/67 computer. For 200 grid points along the 
flexible cylinder length, total run time was approximately 
twelve to fifteen seconds when uSing a precompiled source 
deck. The program is listed after the Two Degree of Freedon 


cascade listing. 
E. RESULTS 


Pressure coefficients and generalized aerodynamic forces 


were computed for sinusoidal deflection modes, 


Z4(X)= LW (x) = sin(j7x) , j = 1,2, ©" (VITI- 35) 
J 


These cases were previously studied by Widnall and 
Dowell, [Ref. 27], who derived a solution of the full 


linearized eguation using a Laplace~transform technique. 


For figures 53 through 55 the grid fineness ( number of 
increments along a Mach line between the inner and outer 
cylinder ) was adjusted as a function of the difference in 
radii (or R/L) to maintain 200 points along the unit 
cylinder length. The inner radius was kept fixed at 0.0002 


while the outer radius was varied fron 0.1 to 5.0. 
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For low values of the circumferential mode number im + 
pressure distributions were excellent, with sharp 
discontinuities evident at reflection points on the outer 
cylinder. However, for n > 2, these reflection points 
exhibited a marked oscillation in pressure whose amplitude 
and duration were dependent on the inner radius and grid 
fineness used. By arbitrarily varying these input 
parameters, the oscillations could be minimized, but not 
entirely removed. Over a small range of inner radii, 0.01 to 


0.0002, overall effect on Q was celative small with 
i 


differences of less than fifteen percent in most cases. For 
R/L large enough that no reflection occurred, Q and 
mo 


pressure distributions were independent of inner radius size 


and grid fineness. used. Further investigation is being 
conducted to develop stability criteria as a function of 


Mach number, grid fineness, and R/L. 





— Widnall and Dowell (Ref. 27) 
o Present study (e n=6) 
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Figure 53. Generalized aerodynamic force Qjj}p versus radius/length 
ratio as a function of the circumferential mode number (Rj=.0002). 
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Figure 54. Generalized aerodynamic force Q}]] versus reduced 
frequency as a function of Mach number (R;=.0002). 
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—— Widnall and Dowell 
(Ref. 27) 


Present study 
M= /2 
= O 


5.0 
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J ‘ 

R/L 
Figure 55. Generalized aerodynamic force Qj2R versus radius/length 
ratio as a function of the circumferential mode number (Rj=.0002). 
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F. CONCLUSIONS 


The method of characteristics solution of the full 
linearized unsteady velocity potential equation has been 
developed to compute the pressure distributions and 
generalized forces for harmonically oscillating cylindrical 
shells with arbitrary radius-to-length ratios, axial and 
circumferential mode numbers, supersonic Mach numbers and 
reduced frequency. In addition, the method is easily 
adaptable to include arbitrary (non-sinusoidal) boundary 
conditions. Moreover, it is felt to be more easily applied 
than the complex integration technigues reguired for the 


Laplace-transform solution. 


V4a4 





The following 
the progran. 


ICO 


ID 


IFLUT 


IHALT 


ITHAVEP 


IPRES1,IPRES2 


ITER 


ITO 


ITRMAX 


APPENDIX A 


TWO DEGREE OF FREEDOM CASCADE 


variables appear in the logic statements of 


used in subroutine GENFPT. When equal to one, 
it indicates that the previous grid point was 
a lower surface (blade or wake) point. ICO is 
always set equal to zero before exiting from 
GENFPT. 

integer array with each element corresponding 
to a cascade blade. Normally zero, elements 
are set equal to one the first time COEF2 is 
called for each blade. 

an input variable. Set equal to one if 
flutter solutions are not desired; zero 
otherwise. 

a termination flag. When equal to one, X, as 
determined from the real Or imaginary 
component of the flutter determinant, is 
negative. 

the number of the computation step along a 
Cight~running Mach line. Set equal to zero at 
the initial left-running Mach line, it is 
incremented until it is equal to LIMIT. 

input variables identifying blades for which 
specific information is desired 

iteration counter for the number of times 
through the flutter routine. 

an information flag. If other than zero on 
return from FLUTER, the discriminant of the 
quadratic for the real roots of X is less 
than zero. 


maximum number of flutter solution iterations 
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IWRITE 


JLINE 


LB 


LCOUNT 


LIMIT 


LINW 


LVEC 


NGRDFN 


NEWDST 


NITOT 


NN 


allowed before program termination. 
output variable. Set equal to one in COEF2 
after integration of the pressure 


distribution for the foe blade is completed. 
the Mach line counter. Incremented each time 
a right-running Mach line is completed, itis 
reset to the value of IHAVEP at the leading 
edge of a new blade. 

an integer array. This is the first subscript 
of the matrix used for storing the discrete 
pressure and incremental moment distributions 
for each blade. 

the value of IHAVEP at a surface or wake grid 
point. When more than one blade is crossed by 
a Cight~running Mach line, LCOUNT is set 
equal to OLDL at the initial left-running 
Mach line and incremented by NGRDFN after 
each blade. 

the Maximum value of IHAVEP along a 
cight-running Mach line. 

the value of LCOUNT at the last grid point on 
a blade. If IHAVEP is greater than LIMW, the 
grid point is in the wake. 

identical to LB, this value is passed to the 
Main program via common. 

an input variable. It is the number of 
increments along a fright~running Mach line 
between blades. 

the number of Ax increments along a blade 
chord (of unit length). 

the total number of grid points on a _ blade. 
This variable is passed to the main program 
via common. ; 
an integer array in COEF used to identify 


successive grid points on a blade. 
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NSTPTS 


NUM 


NUMBLD 
NUMOLD 


NUMPI 
MAXT 


the number of Ax increments between the 
initial Mach line and the leading edge of the 
second blade. 

one less than the blade actually encountered. 
The program numbers blades from zero to 
NUMBED = 7. 

the number of blades in the cascade. 
preserves the maximum value of NUM 
encountered along a right-running Mach line. 
NUMOLD + 1. 

the value of JLINE on a right-running Mach 
line which intersects the leading edge of a 
blade. 


Other program variables. 


AI 
A3,A4 


ANGP12,ANGM12 


ANGP34, ANGH34 
Ba 

CR, CI 

DELTA 

DELTAS 

DR,DI 


DSTSTR 
FAZE 
Beat, Fr 2, FES 


1 
A, the factor —kAy 
I 2 
real and imaginary components of a 
real and imaginary components of Bs 


phase angle between real and imaginary 
components of the lift force and moment due 
to plunge 
Same as for ANGP12,ANGN12, but for pitch 

Me 
Me-1 
CROs « defined by Eqs. (IV~-13c,d) 





1 
B , the factor —k Ax 
L 4 


the interblade phase angle in radians 

the step size along a characteristic (Mach 
line) 

aa defined by Eqs. (IV-13e,f) 


the step size along a blade 
the interblade phase angle in degrees 
the flutter frequency based on *he imaginary 


root and the two real roots respectively 
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Hod,6S2,FS3 


FSTRMN 
Gi, eee G9 
HDSTRL 
RaI2R, ««. 
Teli. eee L& 


K561 


M1, ess M4 


OMEGAA 


OMEGAH 
Pei2,P112 


PR34,P134 
P12MAG, P34HAG 
REDFRQ 

RINV 
RMR12,RMI12_ 
RMR34,RMI34 


R12MAG, R34MAG 


het, RL2 
RL3,RLY 


RM1,RM2 


the flutter speed based on the imaginary root 
and the two real roots respectively 

the input Mach number 

defined 
solves for perturbation quantities 
one-half DSTSTR 


variously in each subroutine which 


defined in IfLI-B for each subroutine 

reali and imaginary lift coefficients for 
plunge and pitch oscillation respectively. 
Defined in Section IV. 

real and imaginary moment coefficients for 


plunge and pitch oscillation respectively. 
Defined in Section IV 
Enowiaetor 7 t= 

w 2 
the factor BG) 
matrix storage for the real and imaginary 
pressure distribution resulting from plunge 
Matrix storage for the real and imaginary 


pressure distribution resulting from pitch 


absolute magnitude of the complex lift force 
due to plunge and pitch rexpectively 

the reduced frequency 

(1/k) 


real and imaginary components of the moment 
distribution due to plunge 

real and imaginary components of the moment 
distribution due to pitch 
of 


resulting from plunge and pitch respectively 


absolute magnitude the complex moment 


real and imaginary component of L 
the real and imaginary components of L 
a 


real and imaginary components of an 
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RM3,RM4 


RMUU 
RSUBA 

S 

STGANG 

STEP 

TNWDST 

TRNGLH 

fe22R, .«.. C221 
U33R, .-. C331 
Mate, ... CHAI 


Ween, «+ COSI 


VRPLNG, VIPLNG 


VRPANL, VIPANL 


XI 
XN 
XNEW 


XLNGTH 
XPT 


real and imaginary components of Me 


the blade density parameter 
the radius of gyration 
the factor /M2-1 
the input stagger angie, in degrees 
amount by which (1/k) is incremented after.a 
call to FLUTER 
input variable defining the vertical distance 
between adjacent blades 
sin (a) eDELTAS. The distance between blades 
divided by the grid fineness ratio 
perturbation quantities resulting from pitch 
at grid points not on a lower surface 
same as above, but at grid points on a lower 
surface 
perturbation quantities resulting from plunge 
at grid points not on a lower surface 
same as above, but for grid points on a lower 
surface 
real and imaginary components of the normal 
velocity perturbation due to plunge at the 
leading edge of the first blade 
real and imaginary components of the normal 
velocity perturbation due to pitch at the 
leading edge of the first blade 
M2 

M24 

the ratio of natural frequencies of plunge to 





the factor 


pitch 

/im {xX} 

floating point value of NUM 

blade position relative to the leading edge 
of the plade 
length of a Mach line between adjacent biades 


matrix array to store the relative blade 
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XR1,XR2 


XSUBA 


XSUBO 


position of each grid point 

the square roots of the two real roots of the 
flutter determinant 

the position of the center of gravity, 
measured from the elastic axis 


the elastic axis position 
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The following 


APPENDIX. B 


PANEL FLUTTER 


Variables appear in the Cylindrical 


ShellysPanel Flutter progran. 


Al 
A2 
A3 


AY 
BETA 


CPR,CPIL 


DELTAS 
DPHIS1 


DPHIS2 


DPHIDX 
DPHIDR 
FACT1 
FACT2 
FSTRAN 


IHAVEP 


IPRINT 


LCOUNT 


AS 
the factor — 
QM 


M 
the factor 


2/N2=-1 
kHAs 


the factor ——— 
2/M2-1 


M 
the ee ey 


fizai 
real and imaginary components of the pressure 
coefficient 

distance along a Mach line between radii 

the directional derivative of g¢ along the 
left-running Mach line 

the directional derivative of g along the 
right-running Mach line 


g 


x 


k2 

n 2 

Sau 

the Mach number 

the number of a computational step along a 
right-running Mach line 

set equal to one if the user desires the 
entire flow field printed out. Otherwise, 
zero 


consecutive number of the grid points on the 
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N 
NGRDFN 


QREAL,QIMAG 


REDFRQ 
RI, RO 
TANR 
XPT 


outer cylinder surface. Incremented by one 
after each call to RAD1 

the axial mode number, am 

mode number of the panel displacement used in 


the computation of Q 
nr 


the circumferential mode number, n 

the number of increments along a 
right-running Mach line between the outer and 
inner cylinder 


real and imaginary components of Q 
mr 


the reduced frequency k 

the inner and outer cylinder radii 

the flow tangency condition 

axial position of grid points on the flexible 
panel surface 


ise Me 
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OOO 


(=) AOQ OOOAQ AAG MAOM AQ OM 


ADO O 


AM 


1 CONTINUE 
INITIALIZE FLOW FIELD QUANTITIES AND LOGIC VARIABLES 
CALL INITAL 


CALL COEF 

FAC = 1./DSTSTR 

NEWDST = FAC 

CIMIT = ( NUMBLD - 1 )*NSTPTS + NEWDST + 1 

LIMW = NEWDS 

MAXI = _NGROFN + NSTPTS 

NUM = 0 

NUMOLD = 0 

NUMPI = 1 

OLDL = 0 

IcQ = 0 

IWRITE = 0 

THE INNER LOOP IS FOR FLOW FIELD CALCULATION 
2 CONTINUE 


CALL COMPXY 


IF POINT IS ON THE INITIAL RIGHT RUNNING MACH LINE, 
CALL MACHLN. 


PEC CCUUNT eeO2 0 5) Go 10 3 
IF THE POINT IS ON A BLADE OR IN A WAKE, GO TO 4 
Tent THAVER 2EQ. £COUNT } GO TO 4 


IF THE POINT IS ON THE INITIAL LEFT RUNNING MACH LINE, 
CALL MACHLN. 


ECORI AVERS.EG. U ) GO TO 3 
CONDITION FOR THE FIRST ENCOUNTER WITH A BLADE 
IF cree eEQ. NUMPI*NSTPTS ) AND. 


1 INE -EQ. MAXI 32 GO TO 10 
CALL GENEFPT 
IF ( THAVEP .EG. LIMIT ) GO TO 7 
IHAVEP = THAVEP + 1 
GO TO 2 
3 CALL MACHLN 
IF ( THAVEP .EQ. LIMIT ) GO TO 7 
IHAVEP = IHAVEP + 1 
GO TO 2 
4 100 = 1 
IF THE POINT [S IN THE WAKE, CALL WAKE 
IF ( IHAVEP .GT. LIMW ) GO TO 5 
CALL TOP 
CALL BOTTOM 
CALL COEF1 
IHAVEP = IHAVEP 1 
IF ( NUM «EQ. O } GOTO 2 
NUM = NUM - 1 
LCOUNT =_LCOUNT + NGRDFN 
LIMW = LIMW — NSTPTS 
GO TO 2 
5 CALL WAKE 
CALL COCF2 
IF IWRITE EQUAL 1 ON RETURN FROM COEF2, IS THE 
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an 


ANAND 


NANO 


ANDAAM AMMO O 


AAO 


QAM 


1 
6 


10 


PRESSURE DISTRIBUTION FOR THIS BLADE DESTRED? 


IF a!  AWRITE oNE. 1 2 GO TO 6 
Ve ( € M «EQ. IPRESL ) .OR-. ( M EQ. IPRES2 ) ) 
GO TO 14 


IWRITE = 0 
IF ( ITHAVEP .EQ. LIMIT ) ICO = O 


IF THIS IS THE LAST_POINT TO BE COMPUTED IN THE FLOW 
FIELD AREA, TERMINATE. ; 


a (( NUM .GE. oe ).AND.CTHAVEP .EQ. LIMIT )) 


GO TO 
IF THIS IS THE LAST GRID POINT ON A RIGHT-RUNNING 
MACH LINE, SWITCH THE_THE LOGIC VARIABLES AND GO TO 
THE NeXT MACH LINE. OTHERWISE: CONTINUE. 
DES CSIHAVEP cEQ. LIMIT 1 GO 10 7 
THAVEP = ITHAVEP + 1 
ITF ( NUM .EQ. O ) GO TO 2 
LCOUNT = LCOQUNT + NGRDFN 
NUM = NUM - 1 
GO TO 2 
CONTINUE 


IS THIS THE INITIAL RIGHT-RUNNING MACH LINE? 
IF ( LCOUNT -EQ. 0 } GO TO 8 


THIS IS_THE LAST GRID POINT ON A RIGHT-RUNNING MACH 
LINE. STORE THE PERTURBATION QUANTITIES JUST COMPUTED 
AND GO TO THE NEXT MACH LINE. 

I = IHAVEP + 1 

U22R(I) = U2RNEW 

U2Z2I(I) = U2INEW 

V22R¢1T) = V2RNEW 

V22T(1) = V2INEW 

C22R(1) = C2RNEW 

C221(1) = C2INEW 

U44R(1) = U4RNEW 

U441I¢1) = U4INEW 

V44R(1) = V4RNEW 

V44I(1) = V4&INEW 

C44R(1) = C4RNEW 

C441(01) = C4INEW 

SWITCH THE LOGIC VARIABLES TO START THE COMPUTATION. 
IHAVEP = 0 

NUM = NUMOLD 

ICO = 

JLINE = JLINE + 1 

IF THIS iS NOT THE FIRST BLADE, RESET THE LOGIC SWITCH 
VARIABLES LCOUNT, LIMW, OLDL. 

IF ( NUMOLD .GT. 0 ) GO TO 9 

LCCUNT = LCOUNT + 1 

GO TO 2 

LCOUNT = OL 

LIMW = oe TE TS + NEWDST 

OWL, = Sot 

GO TG 2 

CALL NEWBLD 

IcO = lL 
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OQ AOOOH 


AO 


aAamM 


ANDO 


ANOO 


Zx=eEeecccr TN SET me 


ADAM -O 


TO = 0 

I = NUMBLD 

ACTGR IS USED TO MAKE LIFT AND MOMENT VALUES AGREE 
ITH GARRICK AND RUBINOW NOTATION. 
ACTOR = 0.5/(REDFRQ*REDFRO) 

L = RLICI1I)*FACTOR 

2 = RL2CII)*FACTOR 

3 = RL3(11)*2.* FACTOR 

4 = RL4&(11)¥*2.*F ACTOR 

L = RMICLIP*2.*FACTOR 

2 = RM2(11)*2.*F ACTOR 

3 = RMZ(L1)*4.*F ACTOR 
M4 = RM4(11)*4.*F ACTOR 


DEFINE CONSTANTS TO BE USED IN THE FLUTER SOLUTION. 
RINV = 1./REDFRQ 


W2 = WHIVA*WHWA 

R2 = RSUBA¥*RSUBA 

X2 = XSUBA*XSUBA 

OMEGAA = RMUU*R2 

OMEGAH = RMUU*W2 

Din= sie = Earls £2245 = LSeMzZ 

DR = Li*¥M3 - L3*M1 - L2*M4 + M2*L4 

Cos—shnUU=( ASUBAS( Mi + £3 3) - € M32 = OMEGAA ) - 

( L1*R2 + RMUU*X2 ) ) + OR 

I = RMUU*( XSUBA*¥( M2 + L4 } - M4 - Le*R2 ) + DI 
F ( WHWA .EQ. 0. ) GO TO 20 
EFINE CONSTANTS FOR QUADRATIC EQUATION OF REAL 
OMPCNENT AND FOR THE IMAGINARY COMPONENT. 


(L1 —- RMUU )/OMEGAH + M3/O0MEGAA - le 
CR/( OMEGAA*OMEGAH) 

L2*OMEGAA + OMEGAH*M4 

SOLVE FOR AND TEST THE DISCRIMINANT 


DISCRM = B*¥B -_4.*C 
IF ( DISCRM eLT. O- 2 GO TO 40 


SOLUTION OF IMAGINARY COMPONENT 
Xl = =2c1/0 
F 

‘ 


I elT.e. O. ) GO TO 41 
RT(XIT) 


I 
Q 
ON OF R 


AL COMPONENT 
QRT(DISCRM) )*0.5 
QRT(DISCRM) )*0.5 
O. ) -OR- ( XR21 
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SIMQ 1S A SIMULTANEOUS LINEAR EQUATION SOLVER 


A - MATRIX (N BY N2) OF COEFFICIENTS STORED 
COLUMNWISE (DESTROYED) 
3 = VECTOR OF GRIGINAL CONSTANTS. FINAL SOLUTION 
ON THIS VECTOR. 
N - NUMBER OF EQUATIONS AND VARIABLES 
MUST BE > 1 
KS - OUTPUT DIGIT 
O FOR NORMAL SOLUTION 
1 FOR A SINGULAR SOLUTION 
DIMENSION A(1),B(02) 
FORWARD SOLUTION 
TOL = 0. 
= 0 
= -N 
65 J = 1,N 
=Jtil 
=JJ+Ne tl 
GA = e 
SB ddo J 
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THIS IS THE MAIN PROGRAM 


COMPLEX PHI(400),DPHIS1(469) ,D0PHIS2(400) 
DIMENSION X(400) ,R(400) »X#T(400),CPR(400),CP1(400) 
DIMENSION DUMMY(400),DUNCE (400) 
COMMON/&6LCKI/ FSTRMN, DELTAS, Dx, DHeDSTSTR yROyRI »REDFRO 
COMMON/ BLCK2/ AlyA2,A3,A%sRNaRMyRMLeFACTL,FACT2 
COMMON/BLCK3/ NGRDFN«NPTS,LCOUNT,THAVEP sNeMeM1,T PRINT 
COMMON/BLCK4/ PHI,CPHIS1,DPHIS2 
GMMGN/BLCK5/ XsR_yXPTeCPR,CPI 
DATA PI/3.141593/ 
Ree uote eDEL A Soeveu= Vio wei XO =AZ )/0ELX 
+ 
EPS = 1.E-04 
TER = 0 
CALL INPUT(IER 
LE (otek 2.60. 1) GO TU 12 
CONTINUE 
THAVEP = 1 
LCOUNT = 1 
CALL INITAL 
CONTINUE 
CALL COMPXR 
IS THE GRID POINT ON THE INITIAL RIGHT RUNNING MACH 
LINE? 
IF (€ CCOUNT .EQ. 1 ) wAND. (CYHAVEP .LE-] NGRDFN )) 
L GO 1G 3 
IS THE GRID POINT ON THE OUTER CYLINDER SURFACE? 
IF ( THAVEP .EQ. 0 ) GO TO 5 
IS THE GRID PGINT ON THE INNER CYLINDER SURFACE? 
IF € ITHAVEP «EQ. NGRDFN ) GO TO 6 
THEN THE GRID POINT IS IN THE GENERAL FLOW FIELD. 
CALL GENFPT(IER) 
IF ( ITER .EQ. } GO TO 12 
IHAVEP = IHAVEP + 1 
GO TO l 
CALL MACHLN 
IF ( THAVEP .EQ. NGRDFN } GO TO 4 
THAVEP = IHAVEP + 1 
6G 10 1 
THIS IS THE LAST POINT ON THE INITIAL MACH LINE. 
LCOUNT = LCOUNT + 1 
IHAVEP = 0 
GO TO 1 


CALL KAZ 
IS THIS GRID POINT PAST THE FLEXIBLE CYLINDER LENGTH? 
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